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Our Heritage from Wilbur and 
Orville Wright 


HUGH L. DRYDEN* 


National Advisory Committee for Aeronautics 


fen FIFTIETH ANNIVERSARY of the historic event at 
Kitty Hawk furnishes the opportunity for in- 
numerable assessinents of the significance of the first 
powered flight to mankind, surveys of the astronomical 
progress of aeronautics during the last 50 years and 
forecast of still greater achievements for the future, 
and simple tributes to the memory of the Wright broth- 
ers. In the name of the aeronautical scientists and en- 
gineers who constitute its authors and readers, the 
JOURNAL OF THE AERONAUTICAL SCIENCES pauses to 
honor the memory of the two men, who, without forinal 
training in science or engineering, were the first success- 
ful aeronautical research scientists and engineers. From 
them we have received a great heritage. 

To many, the greatest gift of the Wrights was the 
airplane itself and all that it has become. On the fourth 
flight on December 17, 1903, Wilbur piloted the first 
airplane for a distance of 852 ft. in 59 see. On the 25th 
anniversary, the official records showed that airplanes 
had remained in the air for 65 hours and 25 min. and 
had traveled 4,466 miles in flight without refueling. 
The record speed was 318 m.p.h. and the record altitude 
38,418 ft. Today, as this note is written, the distance 
record without refueling is 11,236 miles, the record 
speed is expected to be confirmed as 753.4 m.p.h., and 
the official record altitude is 63,668 ft. A research air- 
plane carried aloft in a mother airplane has reached un- 
official records of 1,238 m.p.h. and 83,235 ft. More im- 
portant than these records, the airplane is the key to 
transportation, shrinking the world, and is the dominant 
weapon in the protection of the nation, deterring the 
aggressor. The activity of the Wrights was the well 
spring from which the stream of modern aeronautics 
was generated. 

To aeronautical engineers of the present day, the 
greater heritage is the manner in which the pioneer 


* Director. 


The Wright brothers, Wilbur and Orville. 


aeronautical engineers tackled the problems of their 
time and the qualities of character which contributed to 
their success. The first flight was the accomplishment 
of a team of two men rather than of an individual. Two 
minds are better than one in digesting and interpreting 
the past experience of others and in examining new 
ideas. The importance of the lengthy discussion and 
criticisms of each brother's ideas by the other cannot be 
overemphasized. Although the current mathematical 
theories were not useful, both Orville and Wilbur ap- 
plied their keen intellects to devise physical theories of 
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the behavior of their gliders. The team method of ap- 
proach, the development of ideas in advance of hard- 
ware, and the critical examination of proposals by 
several minds before final decision are as effective today 
in the successful development of new airplanes as they 
were 50 years ago. 

The Wrights were experimenters as well as theorists. 
Their flight experiments with gliders which demon- 
strated errors in the published data, their attempts to 
find suitable laboratory methods, and their successful 
use of a small wind tunnel have been described many 
times. The ingenious manner in which sources of error 
were eliminated or minimized in the wind-tunnel ex- 
periments should be studied by every engineer. Labora- 
tory experiments, however self-consistent, were not re- 
garded by them as conclusive until verified by flight 
measurements. The airplane structure was subjected 
to proof tests. To a surprising degree the design of the 
first airplane was based on a well-rounded synthesis of 
rational thinking, laboratory experiment, and _ flight 
verification. 

The Wrights demonstrated the role of inventive 
ability in technological development. They were in- 
ventors of new methods and devices, but they were 
much more. They were aeronautical engineers and 
designers, stress analysts, practical constructors, and 
businessmen. Inventive ability is a skill whose fruitful- 
ness depends on the environment in which it is exer- 
cised. The professional inventor appears less often 
today in our more complex and highly developed tech- 
nology, because a single individual rarely has an ade- 
quate technical background in more than one field. 
Inventive ability, added to other specialized knowledge 
and skill on the part of the members of the group, is in- 
dispensable in successful development. 

Aside from these lessons in the techniques of our pro- 
fession, we find much food for thought in the spirit and 
motivation of the Wright brothers in their work. In 
his early letters Wilbur characterized his belief that 
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flight was possible to man as a disease that had gripped 
him. He showed clearly that the expectation of fame 
and fortune were not important elements and spoke often 
of the pleasure of scientific investigations. On another 
occasion Orville wrote of ‘‘the fun we got out of learn- 
ing new truths.”’ The urge to blaze new trails, the 
challenge of the unknown, the desire to achieve some 
feat hitherto impossible are still more powerful motiva- 
tions to many than monetary reward. 


The Wrights carried out their fascinating experi- 
ments with modesty, courage, and patience. They did 
not take undue risks and had faith in the ultimate out- 
come of their work. They were, however, completely 
devoted to their work. Though not addicted to gam- 
bling, they did in fact stake their time, their money, their 
reputation, and their very lives on the great adventure 
on which they were engaged. Fortunate is the aero- 
nautical engineer today whose task develops similar 
courage, patience, and faith. 


Wilbur and Orville thought that they were intro- 
ducing an invention that would make wars impossible, 
and they were presented with medals by the French 
Peace Society. Time has taught us otherwise—that 
airplanes are the leading military weapons, as wel! as a 
great social tool in time of peace. Even in this connota- 
tion, however, we hope and work that the desire of the 
Wrights may come true—that the airplane in the hands 
of men of good will may become a tool for the preserva- 
tion of peace. 

The greatest tribute that we can pay to the memory 
of Wilbur and Orville Wright is further to advance the 
cause to which they devoted their own lives, to push 
back still further the barriers of the unknown, to con- 
vert our dreams of progress to accomplishment, im- 
pelled forward by what the human mind can conceive. 
Such a memorial is more lasting than the granite that 
lies in the sands of Kitty Hawk to mark the site of the 
first powered flight. 
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The Compressible Laminar Boundary Laver 
with Arbitrary Pressure and Surface 
Temperature Gradients’ 


DEANE N. MORRIS? ann JOHN W. SMITH? 
Bell Aircraft Corporation 


ABSTRACT 


The Karman-Pohlhausen method is extended systematically to 
cover the hitherto unsolved problem of two-dimensional com- 
pressible laminar flow with arbitrary surface temperature and 
pressure gradients and arbitrary but constant Prandtl] Number. 
The main features of this development are clearly distinguished 
from special numerical assumptions (e.g., particular boundary- 
layer profile representations), which are evaluated a posteriori and 
are subject to modification and refinement. Application of the 
present method to the calculation of heat transfer and skin fric- 
tion in a given case involves the numerical integration of two 
simultaneous ordinary differential equations, which is readily 
accomplished by hand or with an automatic digital calculator. 
The theory is well adapted to an investigation of the combined 
effect of pressure and surface temperature gradients in the com- 
pressible boundary layer, for which no other practical method is 
available. 

Solutions by the present method are compared with more exact 
solutions in the special cases for which the latter are available, and 
the agreement is good in all cases. The special cases used for 
comparison include skin friction for incompressible retarded flow, 
skin friction and heat transfer for incompressible flow over a spe- 
cial cylindrical shape, and skin friction and heat transfer for 
compressible flow over a flat plate, with and without a surface 
temperature gradient. 

A sample problem is presented for a circular-are airfoil at zero 
angle of attack in supersonic flow, with a specified surface tem- 
perature distribution. The local skin-friction and heat-transfer 
coefficients are determined as a function of distance along the sur- 
face and are compared with the results obtained by application of 
the common flat-plate relations. 


SYMBOLS 


a = local speed of sound 

oy = local skin-friction coefficient 

Cp = specific heat at constant pressure 

f(%,) = velocity profile 

fo fi = functions of y occurring in velocity profile 
= temperature profile 

%, £1, 2 = functions of yn occurring in temperature profile 
hy = local heat-transfer coefficient 

HA*) = universal functions of A* 

k = thermal! conductivity 

l = constant characteristic length 

M = Mach Number 


Received July 28, 1953. 

* This work was done under Air Force Contract No. W33- 
038-ac20063, Project MX-1677B, sponsored by the Equipment 
Laboratory, Wright Air Development Center. 

t Project Aerodynamicist. 

t Aerodynamicist. 


805 


Nuz 
Nuc 

Pr 
qu 


local Nusselt Number 

Nusselt Number = 
pressure 

Prandtl Number 

local rate of heat flow at the wal 
heat flow parameter 

recovery factor 

airfoil radius of curvature 


gas constant 

Reynolds Number = 
local Reynolds Number = tox/v 
free-stream Reynolds Number = 

distance measured from center of airfoil 
Sutherland constant = 198°R 

absolute temperature 

velocity components in the x, y, 2 directions, re- 


spectively 

1/[1 + — 1)Mo?] 

velocity parameter defined in Eq. (11) 

viscosity parameters [see Eqs. (13) and (18) ] 

curvilinear coordinate along the body surface in the 
local free-stream direction, with origin at the 
leading edge 

curvilinear coordinate perpendicular to the body 
surface, directed into the stream 

curvilinear coordinate along body surface per- 
pendicular to the local free-stream direction 

constants, defined by Eq. (26) 

ratio of specific heats 

value of at y = 

dynamic boundary-layer thickness 

thickness parameter defined by Eq. (8 

value of yn at y = Ay 

thermal boundary-layer thickness 

thickness parameter defined by Eq. (8) 

A/é 

slope of surface of airfoil 

variable defined by Eq. (5) 

n/A 

thickness parameter defined by Eq. (8) 

thickness parameter defined by Eq. (8) 

velocity profile form factors defined by Eqs. (13) 
and (14) 

Te" =—1 

temperature profile form factors defined by Eq. 
(18) 

absolute viscosity 

kinematic viscosity = p/p 

density 

shear stress at the body surface 


= 
= 
| R = 
Ro = 
| ke = 
| = 
PY = 
= 
S* = 
T = 
| 
| “vw = 
| Ny? = 
Va- | 
= 
X1, X2 
Day x = 
the 
ish | 
= 
m- 
ve. | 
at | = 
he 7 = 
= 
| 5, = 
| 
A = 
Ay = 
= 
= 
€ = 
=> 
n = 
= 
0 = 
ML, = 
Ato = 
| Aa, Are = 
= 
= 
p = 
Tu = 
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Sub- and Superscripts 


om = quantity made dimensionless through division by a 
constant—e.g., velocity by ~/2ia:, distance by 
/, temperature by 7'~;, density by p.,, etc. 

a = quantity made dimensionless through division by 
function of x—e.g., velocity by mo, distance by 
6, temperature by 7), density by po, etc. 


ce = differentiation with respect to 7 or 7 
i = differentiation with respect to x or & 
t = stagnation conditions 
= free-stream conditions 
ot = free-stream stagnation conditions 
0 = local free-stream conditions (i.e., outer edge of 
boundary layer) 
w = conditions at body surface 
aw = adiabatic wall conditions (insulated surface) 
INTRODUCTION 


PRANDTL BOUNDARY-LAYER EQUATIONS, which 
are derived from the equations of continuity, mo- 
tion, and energy for viscous fluids on the basis of certain 
assumptions,' form the starting point for almost all of 
the present theories of the laminar boundary layer. 
Certain exact solutions of these equations, exact in the 
sense that no further mathematical assumptions are 
imposed, exist in the literature. For incompressible 
flow, the Blasius solution remains as the classic analysis 
for skin friction on a flat plate,' while other exact 
solutions are available for flows with various specific 
pressure gradients,' and for the corresponding thermal 
cases.'~* Exact compressible flow solutions for skin 
friction and/or heat transfer have been presented by 
Crocco,® © Howarth,’ Cope and Hartree,’ and Chap- 
man and Rubesin.® 

All of these exact methods suffer because of the severe 
physical restrictions imposed on the flow conditions 
(e.g., a specific or zero pressure gradient, no heat trans- 
fer, linear viscosity law, etc.), and much effort has been 
expended in the development of approximate solutions 
of the Prandtl equations which will be applicable to a 
more general class of problems." ~~" In most cases, 
these approximations make use of the Karméan-Pohl- 
hausen integral method or some modification thereof. 
In this method, the partial differential equations of 
Prandtl are first integrated across the thickness of the 
boundary layer, and the resulting ordinary differential 
equations are then solved by assuming some particular 
form for the velocity profile. Obviously, these solu- 
tions will not satisfy the original equations except in 
the average across the boundary layer, but by imposing 
sufficient boundary conditions on the resultant pro- 
files, acceptable accuracy can usually be obtained. 
Libby, Morduchow, and Bloom have presented an 
excellent critique of the Karman-Pohlhausen method 
and its modifications. 

In nearly all of the approximate methods cited above 
it is assumed either that the flow is incompressible or 
that there is no heat transfer and the Prandtl Number 
is unity. Either of these assumptions allows the sep- 


aration of the momentum and energy equations, re- 
sulting in a great simplification of the solution, but these 
assumptions also restrict the generality of the results, 
The purpose of this paper is to present an approximate 
but general solution of the Prandtl boundary-layer 
equations, by developing the basic concept of the K4r- 
man-Pohlhausen method in a more systematic manner 
than has hitherto been accomplished. The salient fea- 
ture of this development lies in the fact that the inte- 
grated momentum and energy equations are solved 
simultaneously, without the aid of further assumptions 
or restrictions. The resultant solutions give a clear 
indication of the effects of pressure gradient, surface 
temperature gradient, and their interactions on the com- 
pressible laminar boundary layer. 

Emphasis has been placed on deriving a set of equa- 
tions which is readily solvable by a numerical process 
rather than on obtaining explicit solutions. A clear dis- 
tinction between the essential theory and special as- 
sumptions of a numerical nature has been preserved 
throughout; all assumptions that are basic to the theory 
have been placed in brackets and designated by Roman 
numerals for easy identification. 


DEVELOPMENT OF ESSENTIAL EQUATIONS 


The Prandtl boundary-layer equations for compres- 
sible flow are generally expressed in the form 


(0Ox) (pu) + (0/dy) (pv) = 0 (1) 

Ou Ou op 4 ( 9) 

(2 
Ox Ov Ov Ov Ov 

dp dy = 0 (3) 


u v = + 
ox Ov C ov + 
*) (4) 
c, \Ov 
These are derived from the equations of continuity, 


motion, and energy for viscous fluids,' with the aid of 
the following assumptions: 


[Steady flow (0/0/ = 0)] (I) 
[Negligible body forces (F, = F, = F. = 0)] UD 
{Two-dimensional flow (0/0z = 0)] (IIT) 
[Prandtl boundary-layer assumptions | (IV) 
[Negligible curvature of surface | (V) 
[Perfect gas (pb = pRT)] (VI) 
[Constant specific heat (c, = constant) ] (VII) 


The “Prandtl boundary-layer assumptions” (IV) in- 
volve assuming that all viscous effects are concentrated 
in a thin layer next to the surface, that higher orders 
of the nominal thickness of this layer are negligible, and 
that the Reynolds Number of flow is large. The latter 
not only implies that the viscosity of the fluid is small 
but that the extent of the boundary layer along the sur- 
face is large with respect to its thickness. 
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The so-called integrated momentum and energy 
equations are derived by integrating Eqs. (2) and (4) 
across the boundary layer. Although exact solutions 
of the boundary-layer equations indicate that the ve- 
locity and temperature in the boundary layer approach 
the free-stream conditions asymptotically,! finite thick- 
nesses can be defined for these layers by specifying the 
distance from the surface at which the relative defect 
is less than some given value. Without indicating the 
actual magnitude of this limit, it will be assumed that 


thiscan be done. Thus there exist 


[Finite dynamic and thermal boundary-layer thick- 


nesses, 6, and A,] (VIII) 


The boundary conditions for the dynamic and thermal 

boundary layers are 

= 0, u=v=0 
= Ou/ov = 0, 


v u= 
y>A, Tt=1 


y 4, 


It is known that the integrated equations assume a 
particularly simple form upon introducing the Dorod- 
nitsvn variable’ y, which is defined by the transforma- 


tion 
dn = (p dy) 
n(0) = 0 
(2) 
n(6,) = 6 
n(A,) = A 


In the integration of the energy equation, further sim- 
plification results from the use of stagnation temper- 
ature as the dependent variable, and the following as- 
sumption is made: 


[Pr = c,u/k is constant] (IX) 


In view of assumption (VII), this is equivalent to as- 
suming that the thermal conductivity is directly pro- 
portional to the absolute viscosity. 

The resulting integrated equations are”’ 


dd Uy’ Tu 
(6, + A, + 29) = (6) 
dx Uy Po 
dd ty’ du 
4+ — 


where 
= [u(Ou Oy) ],, 


= = u*(1 — u*) dn | 
0 
f (1 — u*) dn 
0 
A 
A, (T* — 1) dyn 
J0 
6 = 60* -/ u*(1,— 1) dn 
0 


— [R(OT /dy) 


= 
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In order to integrate Eqs. (6) and (7), it is necessary to 
reduce the six dependent quantities of Eq. (8) to a set 
of two dependent variables. The Pohlhausen method 
accomplishes this by choosing appropriate representa- 
tions for the boundary-layer profiles, by means of which 
the six quantities in question are expressed explicitly 
in terms of two variables—essentially 6 and A. It may 
be pointed out that the profile representations are sig- 
nificant only insofar as they relate the quantities defined 
by Eq. (8). 

The obvious advantage of the Pohlhausen method is 
that the solution of the partial differential equations of 
the boundary layer is replaced by the solution of a set 
of simultaneous ordinary differential equations, which 
in many cases are independent. Because of this, the 
method has been used extensively and, with proper ap- 
plication, can give surprisingly good results. 


Representation of the Velocity Profile 


Let the velocity profile u* = 1/u be represented by 
a function f(x, 7), where 7 = 7/6. Then f(x, 4) satisfies 
boundary conditions of the form 


f(x, 0) = 0 
f(x, 1) = 1 
f(x, 1) = 0 (9) 


f(x, 1) = 0 
where the exact number of conditions at the outer edge 
of the boundary layer is somewhat arbitrary. It will 
be assumed that f(x, 7) can be written in the form 


[u* = f(x, 4) = fola) + | (X) 
where 
fo(0) = 0, fi(0) = 0 
70) = 0, = -1 
fo(1) = 1, fil) = 0 (10) 
fo(1) a 0, A) = 
fol) = 0, fi(1) =0 


Strictly speaking, an infinite number of form parameters 
\,(x) is required to represent the exact solution of the 
Prandtl boundary-layer equations, provided the latter 
exists. The accuracy of the present method can be 
increased by adding one or more form parameters in the 
representation of the boundary-layer profiles. 

Particular functions fo(7) and f;(7) will have to be 
assumed in order to derive numerical results from the 
present theory. Questions relating to the choice of 
these functions will be discussed in the third part of this 
paper, after all of the fundamental equations have been 
developed. 

The relation between \,(xv) and the quantities 6 and 
A can be derived from the boundary-layer equations. 
For this purpose Eq. (2) is written in terms of the vari- 
It becomes 


ables (x, 7). 
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Ox On 
On On 
where 


= pd/por + u(On/dx), 
According to assumption (X) and Eq. (10), 


M(x) = —F(x, 0) = (12) 

Uo On? n=0 

Eqs. (11) and (12) can be solved for \,(x), giving 
M(x) = Aw + Xi f(x, (x, 0)/A* (13) 


where 


or") 
(u*/T*) 


14 


= 
The temperature term g(x, 0)/A* is expressed in the 
notation of the following section. 


Representation of the Temperature Profile 


The temperature profile (7* — 1) can be handled in 
a similar way, except that a greater number of form 


parameters is required. It will be assumed that 


[((7* — 1) = glx, 4.) = + 


where 4, = n/A, Aw = T,* — 1, and 


g(0) = 1, gi(0) = 0, g(0) = O 
#(0) = 0, ¢,(0) = —1, = 0 
(0) = 0, ¢,(0) = 0, = 1 
g(1) = 0, = 0, g(1) = 0) (15) 
&(1) = 0, (1) = 0, (1) = 0 
(1) = Q, @,(1) = 0, go(1) = 0) 


As will be seen presently, the form parameter \ (x) is re- 
quired to account for the effect of surface temperature 
gradients on the thermal boundary layer. 

Both Ap(v) and A~(x) can be related to the quantities 
6 and A by use of the Prandtl equation for the thermal 
boundary layer. In the (x, )-plane, Eq. (4) becomes 


Or, Poll — 
Ox on Pr 
(1— Pr) (16 
dnl TL dn c, On 


According to assumption (XI) and Eq. (15), 


OT | 
An(x) = — g(x, 0) = 


3 BT 
To n=0 


By means of Eqs. (16) and (17), the desired expression 
for \,(x) can be derived directly. To obtain a solution 
for \p(x), however, it is necessary to differentiate Eq. 
(16) once with respect to 7. There result the expres- 
sions 


(17) 
A2(x) 


tig” 


ha = An*A*? = [f(#, 0) + 
1- 


22,9) Fi... 
xf A* ] 


| Tot 

| 0) + SX\Pr — 

1 — — 

0) 
A* | 


where A* = A/6 and 


‘ity (u*/T*) w 


General Solution 


In the preceding paragraphs the form parameters 
of the dynamic and thermal boundary-layer profiles 
have been related implicitly to the quantities 6 and A. 
In this section the six quantities 7,,, 3, 6, and 
occurring in Eqs. (6) and (7) will be expressed in terms 
of Ay, A,,, and A*. This makes it possible to solve 
Eqs. (6) and (7) by a process of numerical integra- 
tion. 

It is convenient for calculative purposes to reduce all 
equations to dimensionless form. For a characteristic 
length, /, a Reynolds Number is defined as Ry = uol'v<». 
A bar denotes lengths made dimensionless by dividing 
by /, and a prime on a dimensionless quantity is under- 
stood to indicate differentiation with respect to 2. 

The quantities 7, and gq, are eliminated forthwith 
from the differential equations by means of Eq. (S), 
and the result is rendered dimensionless. The inte- 
grated momentum and energy equations assume the fol- 


lowing final form :”° 


d ‘ /(y-)) to’ 6,* A,* 3 7 | 
di (Rod*) =2 ty) 0) + + 2 (Rod?) (19) 
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1 2 
(Ry?) 
dx 


where the asterisks on J, 6,, A;, and @ denote division 
by 6. These four quantities can be evaluated by carry- 
ing out the integrations indicated in Eq. (8), using the 
given profile representations. There result the equa- 


tions 
= Bo + Biri + BoA,” (21) 
= Bs + (22) 
Ai* = + + (23) 


A* = A*) An*(A*) + 
+ Ar [Aws73(A*) + 
+ + 
+ + It (24) 
where 
| 4 
Ho( A*) = j :) dj + 
0 A 
A* 7 
H,(A*) = af fo(a) = 


H(A*) = ) dn + (25) 
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A* 7 
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1 
(Rie? (20) 


“1 
= £0 dit 
| 
= gi diy 
° 
0 
= -f — fo?) da (26 Cont.) 
| 
By = f — 3fo?) da 
0 
Bro = af dy | 
| 


1 
0 


For given functions f,(4) and g,(4,), the quantities 
8, and /7,(A*) can be determined. 


Initial Conditions 


In order to integrate the system of differential equa- 
tions which has been derived, it is necessary to estab- 
lish the boundary conditions at the leading edge of the 
surface. In general, this determines initial values for 
\, and A* which will be functions of the velocity and 
temperature at the leading edge. 

According to assumption (IV), the boundary-layer 
equations are not valid near the leading edge of a body, 
because the assumption of large Reynolds Number is 
violated in that region. Two types of flow can exist 
at the leading edge: that for which the velocity is 
zero—i.e., a Stagnation point—and that for which the 
velocity is finite, as for a flat plate or a body in pure 
supersonic flow. In the case of stagnation point flow, 
the exact solution of the Prandtl boundary-layer equa- 
tions also satisfies the complete Navier-Stokes equations 
when the flow is incompressible and has constant vis- 
cosity.' In addition, Illingworth?! has solved the com- 
pressible boundary-layer equations for stagnation-point 
flow, and his results agree reasonably well with the in- 
compressible solution. It can be concluded that, de- 
spite the violation of one of the basic assumptions, the 
boundary-layer equations for compressible flow may be 
applied at a stagnation point. The conditions near 
the leading edge of a flat plate have been investigated 
theoretically by several authors.**~*4 Although these 
more exact solutions differ from the boundary-layer 
solution very near the leading edge, the results of these 
studies indicate that the leading-edge solution ap- 
proaches the boundary-layer solutions asymptotically 
and that at an extremely short distance from the lead- 
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dn 
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ing edge (R, ~ 100) the two are almost indistinguish- 
able. Therefore, the solution of the Prandtl boundary- 
layer equations corresponding to the ordinary initial 
conditions is valid everywhere except in the immediate 
neighborhood of the leading edge. — 

For stagnation-point flow, Ry vanishes at the leading 
edge, and the quantities Rod? and R,6? are zero. For 
flow with finite leading-edge velocity, it will be assumed 
that the dynamic and thermal boundary layers orig- 
inate at the leading edge of the surface—i.e., 


[6 = A = Oat & = 0 for finite edge velocity] (XII) 


As a result of this assumption, the quantities Rod? and 
R,f? are zero at the leading edge for finite edge velocity 
also. Therefore, the following representations are 
valid for any type of flow in the neighborhood of the 
leading edge: 


Rod? = [(d/dz) (Rod) 

Rob? = [(d/d®) (Rob) 
In the limit, the combination of these expressions leads 
to the equation 


d 
— (Rb?) — 0*? (Rod) =(0 (28) 
dz dt z=0 
Eliminating the velocity gradient from Eqs. (19) and 
(20) by use of Eq. (14) and substituting the results into 


Eq. (28), one obtains 


g(t, 0) 


Ato 


(6,:* + + « 
* 

1 — 


+ 0) 
(29) 


This expression serves as one of the two equations 
that determine the initial values of the parameters ), 
and A* for any type of flow at the leading edge. The 
second equation must be obtained from a consideration 
of the specific conditions for finite edge velocity and 
stagnation-point flow. 

(1) Finite Velocity at the Leading Edge. 
one equation has been derived [Eq. (29)] which can be 
used in the determination of the initial values of \; and 
A* for the case where the velocity at the leading edge is 
finite. The second equation that is required for this 
purpose is obtained by combining assumption (XII) 
with Eq. (14) to give 


So far only 


(M10) =9 = 0 (30) 
Eq. (29) now reduces to the form 
— : = 0 31) 
Ex + 1 No” 


These two conditions are sufficient to determine \; and 
A* at the leading edge as functions of Pr, fo, and 7,,*. 

It is apparent that Eq. (31) has a singularity at 1/) = 
0 and T,,* = 1, since £(%, 0) and 6* both vanish identi- 
callv in this case. The condition that the solutions \, 
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and A* exist and be continuous at the singularity 
leads to an algebraic eigenvalue problem, as will be 
shown. 

Neglecting terms of the second and higher orders in 
(T.,* — 1) and 2, one derives from Eqs. (30) and (31 
the following relation, valid in the neighborhood of 
My = Gand T,* = 1: 


( ) | + fo(0)Ho(A*) 
tin? [2f02(0) By (0) + 
Bsf(0) = (32 


The condition that A* exist and be continuous at the 
singularity implies that the equations 


| = (33 
&=0 


Pr fo(0) 


= 0 


=() 


[ Bs 
JS 


be compatible. This certainly will not be the case for 
arbitrary boundary-layer profile functions. 

Suppose gi(7,, 4) represents a family of profile func- 
tions gi(7,), each satisfying the boundary conditions of 
Eq. (15), such that %,(0, \) = A. This generates a 


family of functions //,;(A*, \). It appears now that \ 
functions as an eigenvalue in Eq. (31), arising from the 
condition that A* exist everywhere and be continuous. 
In order to determine }, let the solution of Eq. (35) | 
be represented by Ao*. Substituting this value of 


A* 11 


1 
This 
func 
the 

cond 
and 
pera 
profi 
(2 
tion 
thicl 
and 
resu 


ld 


E 
poin 
(38) 


is an 


= 
} 
| 
a» 
| 
ava 
~ 
1s 
I. 
qd 
1 
I 
\2P 
| 
| | | 


4P6 


arity 


ll be 


rs in 
(31 
d of 


COMPRESSIBLE 


into Eq. (34), there results 


| | Bs 2Bor 
H (Ay d) + + = 0) 
This is an algebraic equation that determines \ as a 


The solution is shown in Fig. 1, using 
Under these 


(35) 


function of Pr. 
the profile functions discussed later. 
conditions, the initial value of A* is given by Eg. (33) 
and is independent of Mach Number and wall tem- 
perature; this result is plotted in Fig. 2 for the same 
profile functions. 

(2) Stagnation Point. 
tion of the boundary-layer equations, the momentum 
thickness J has a finite value at a stagnation point, 
and its derivative is finite. In accordance with these 
results, it will be assumed that 


According to the exact solu- 


(di dz is finite at a stagnation point] (XIII) 


The integrated momentum equation, Eq. (6), gives 
a value for this derivative which, in dimensionless form, 
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+ A* + 20") (36) 


At a stagnation point, Ry = 0, and in order that the 

derivative be finite, the quantity in brackets must van- 
ish, 

+ Ai* + — 0) = 0 (37) 

This expression may be substituted into Eq. (29), 

which reduces the latter to the simple form 

Ty * 0) 

* ’ 90 

[avo + Pr At - = 0 (38) 

For any given Prandtl Number and wall temper- 

ature ratio, Eqs. (37) and (38) can be solved simul- 

taneously to give \; and A* at the leading edge. The 

results are shown in Fig. 3 for the designated profile 


functions. 


Eq. (38) is structurally similar to Eq. (31), and one encounters the same difficulties in the neighborhood of the 
The problem is again solved by treating the parameter {(0) as an eigenvalue of Eq. 


(Ao*) + + 4P [He(Ao*) + ] + 
2/ (8, 0) Ao* [f(%, 0):(0) + 2A, + (8s te + + (39) 


is analogous to Eq. (35) and can be used to calculate £2(0) as a function of Pr, as indicated in Fig. 1. 
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Skin-Friction and Heat-Transfer Relations 


It is customary to express the skin-friction and heat- 
transfer characteristics of the boundary layer in terms 
of the skin-friction coefficient and Nusselt Number, 
respectively. These dimensionless parameters are re- 
lated to the fundamental quantities 7, and g, through 
the equations 


Cy = pollo” (40) 
Nutz = — Tay) (41) 


Here 7 \,,, refers to the adiabatic wall temperature, 
which is defined as the equilibrium temperature as- 
sumed by an insulated surface. The concept of Nus- 
selt Number is useful when the local rate of heat flow 
depends on the local boundary conditions only. In that 
case the adiabatic wall temperature can be calculated 
readily as a function of the local Mach Number and 
Prandtl Number, and the quantity Nu,/1/R, turns 
out to be a slowly varying function of Mo, Pr, and 7,,*. 


This condition is satisfied only on a flat plate with con- . 


stant wall temperature, because this case is distinguished 
by the fact that the boundary-layer profiles are all simi- 
lar—i.e., depend only on the variable y/+/x. When 
the aforesaid condition is not satisfied, then there can 
exist points along the surface where the Nusselt Num- 
ber assumes both positive and negative infinite values. 
This is due to the fact that the temperature 7), at which 
Go vanishes locally does not generally coincide with 
the adiabatic wall temperature, which was derived on 
the assumption that g, vanishes everywhere. To 
avoid this difficulty, it may be necessary to replace the 
Nusselt Number by a new dimensionless heat-transfer 
parameter. For lack of a better choice, the quantity 


= Roly = = (42) 


can be used. 
The given skin-friction and heat-transfer parameters 
can be evaluated by means of the following equations: 


= 2 (5) | (48) 


A*/(#, 0) 


(44) 


SPECIAL FUNCTIONS 


Viscosity-Temperature Relation 


In order to solve the equations of the preceding sec- 
tion, it is necessary to specify a functional relation be- 
tween the coefficient of viscosity and the temperature. 
Of the various viscosity formulas that are currently 
considered in boundary-layer theory, the Sutherland 
law is known to represent the empirical data with the 
greatest accuracy. It has been applied in the present 


theory to the calculation of skin friction and heat 
transfer on a flat plate in compressible flow. The skin- 
friction parameters calculated in this manner agreed 


well with the calculations of Crocco,*®:* but the results 
for heat transfer were not satisfactory. This difficulty 
can be traced to terms involving the quantities Y, and 
X», which become excessively large with increasing rate 
of heat transfer (i.e., large g(%, 0)). The causes of this 
nonconvergence are not fully understood at this time, 
but it is known that the difficulty disappears when the 
above-mentioned terms are suppressed. This is actu- 
ally equivalent to the use of the linear viscosity law, 
u* = CT*, which has been assumed by Chapman and 
Rubesin® and other writers. According to Eqs. (13) 
and (18), this assumption becomes increasingly justi- 
fied as adiabatic wall conditions are approached, be- 
cause the terms in question vanish identically when 
2(%, 0) = 0. This explains the fact that the linear 
viscosity law gives best results for quasi-adiabatic 
surface conditions. 

In view of this discussion, the Chapman-Rubesin 
linear viscosity law will be adopted for the remainder of 


this paper, 
1+S* 


where S* = S’7) and S = Sutherland constant. 


Boundary-Layer Profile Function 


The choice of the functions to represent the velocity 
and temperature profiles depends entirely on the agree- 
ment of the results with the exact solutions. In keep- 
ing with standard procedure for the Pohlhausen-type 
solution, the functions f; and g,; will be chosen as poly- 
nomials. 

Fourth-, sixth-, and eighth-order polynomials have 
been used to represent the velocity profile functions f, 
and it was found that each representation is best suited 
for certain free-stream velocity distributions. For 
monotonically retarded flow, sixth-order or higher 
order polynomials are necessary in the velocity profile 
to assure reasonable accuracy near the separation point, 
and the best results were obtained with eighth-order 
polynomials. For the flat plate, the greatest accuracy 
is achieved using sixth-order polynomials, while at a 
stagnation point or in a region of large acceleration it 
has been found that polynomials of order higher than 
the fourth do not give a solution for the velocity pro- 
file.* In the case of the temperature profile function, 
satisfactory results have been achieved by use of sixth- 
order polynomials. 


It should be emphasized, however, that these results. 


are not necessarily indicative of the maximum accuracy 
obtainable from the present theory. The problem of 
determining the optimum profile functions remains to 
be investigated. 

According to the previous discussion of initial condi- 
tions, the functions g;(7,) and go(7,) must satisfy a spe- 
cial condition in order that the solutions of Eqs. (31) 
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and (38) exist everywhere and be continuous. The = galt) + £0) = 1,2) (46) 
problem is conveniently solved by representing the 
functions in question in the form where g;,(7,) satisfy the conditions 
gu(0) 0, £12(0) = 0, = 0, £22(0) = 0 
£u(0) = 0, £2(0) = 1, g1(0) = 0, g2(0) = 1 
gu(0) = —1, g(0) = 0, g(0) = 0, g(0) = 
= 0, £20) = 0, = 1, g»(0) = 0 (47) 
gu(1) = 0, g2(1) = 0, gu(1) = 0, gee(1) = 0 
gu(1) = 0, f2(1) = 0, gn(1) = 0, go(1) = 
gu(1) = 0, £12(1) = 0, &(1) = 0, f22(1) = 0 
Now & (0) and £(0) can be varied independently, and = 24 — 149° + 28q° — 2097 + 
the functions //,(A*), [/,(A*), and /7;(A*) as-  fi(q) = (1/14) (89 — 74? + 354° — 639° + (51) 


sume the form 


Consequently, Eqs. (35) and (39) reduce to algebraic 
equations that can be solved for the quantities ¢)(0) 
and %(0), respectively. These turn out to be slowly 
varying functions of Pr, as shown in Fig. 1. 

The various methods, based on the use of different 
polynomials, will be designated by the following nota- 
tion: a P will indicate the use of the Pohlhausen method; 
the number preceding the P will indicate the type of 
velocity profile; the number following the P will indi- 
cate the type of temperature profile. Thus, “NPM” 
indicates that the Nth order polynomials were used in 
the velocity profiles and Mth order polynomials in the 
temperature profiles. 

It may be pointed out that the quantity g(0) must 
be calculated by means of a 4P method. This is due 
to the well-known fact that Eq. (37) does not yield a 
solution for \; when higher order polynomials are used 
in the velocity profiles. 

All profile functions satisfy the basic boundary con- 
ditions set down in Eqs. (10) and (15), while the poly- 
nomials g;,(4,) are defined by Eq. (47). Additional 
boundary conditions must be imposed on the third 
derivatives of any profile function to obtain a poly- 
nomial of higher order. Eqs. (49), (50), and (51) list 
the fourth-, sixth-, and eighth-order polynomials, re- 


(48) 


spectively, 

fola) = 24 — 24° + 

fila) = (1/6) — 3%? + — 
fila) = 24 — + 6H — 298 
fila) = (1/10) (24 — 54? + 1094 

104° + 37°) 

= 1 — 2, + — CHF + 208 (50) 
gulf) = 1 + 6G — Sn? + 


= (1/6) (42 — 344 + 392 — 2°) 
Nt 


— 104.4 + 159° — 


4277 — 107°) 


Using the fourth- and sixth-order polynomials, the 
universal constants and functions 6; and //; of Eqs. 
(25) and (26) have been calculated. 


COMPARISON WITH EXAcT THEORIES 


To establish the validity of the special assumptions 
inherent in the method presented above, it was applied 
to the solution of a few problems for which exact solu- 
tions are available. The term ‘‘exact’’ is used here in a 
particular sense—i.e., to designate mathematically ex- 
act solutions of the Prandtl boundary-layer equations 
[(1)-(4)] satisfying some given equation of state and vis- 
cosity law. Numerical solutions in which the error can 
be controlled will also be included in this class. The 
question as to how accurately the physical flow condi- 
tions are described by the boundary-layer equations or 
by the given equation of state or viscosity law consti- 
tutes a separate problem. 

The numerical comparisons referred to above fall into 
two main groups: incompressible flow with positive and 
negative external pressure gradients, and compressible 
flow on a flat plate. 


Incompressible Flow 


The word “incompressible’’ usually implies that the 
Mach Number is zero—i.e., either that the velocity is 
zero or that the speed of sound is infinite. The first 
alternative must be ruled out, of course, which leaves 
the second condition to be satisfied. This cannot be 
done for a perfect gas, since the speed of sound is then 
given by the relation 

a=V yRT 


where y, R, and JT are necessarily finite. It appears, 
therefore, that a perfect gas cannot be considered in- 
compressible in the strict sense of the word and that 
one or the other of these conditions must be relaxed, if 


only slightly. 
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Fic. 5. 


A more accurate definition of incompressible flow, 
with regard to aerodynamics of the boundary layer, 
would state that it is the flow of a perfect gas for which 
the Mach Number effects are negligible. This then 
follows logically as the limiting case of the general com- 
pressible flow problem. 


DECEMBER, 


1958 


However, the type of ‘incompressible’ flow usually 
assumed in boundary-layer theory may be called “‘con- 
stant-density”’ flow to distinguish it from the “low- 
speed”’ flow defined above. For this flow it is assumed 
that the density (and viscosity) are constants, inde- 
pendent of pressure and temperature. The compres- 
sible flow equations constituting the present theory 
[Eqs. (19)-(24)] do not reduce to the correct equations 
describing constant density flow when the Mach Num- 
ber is set equal to zero, because the equation of state 
for a perfect gas has been assumed at certain steps of 
the derivation. This can be corrected by using the 
constant-density flow relations p* = — = 1 wherever 
necessary. 

With this modification, the present theory has been 
applied to the problems of Howarth” (retarded flow 
and Goland‘ (accelerated flow), for which exact con- 
stant-density flow solutions exist. The free-stream 
velocity distributions are 

ty = by — 
Mo = — 
respectively. In the first case only a solution for skin 
friction is available, whereas exact solutions for both 
skin friction and heat transfer exist in the second case. 
The numerical results obtained from the present theory 
are compared with these exact solutions in Figs. 4-6. 


Compressible Flow on a Flat Plate 


In view of the use of Eq. (45), it is proper to com- 
pare the present theory with the theory of Chapman 


9 EXACT SOLUTION 
y 


PRESENT METHOD (4P) 


Fic. 6. Heat-transfer parameter for Goland 
problem (constant density flow, Pr = 1). 
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and Rubesin for the case of compressible flow on a 
flat plate, since both are based on the same viscosity 
law. 

For flows with constant surface temperature, the 
theory of Chapman and Rubesin leads to the following 
results, so far as the skin-friction and heat-transfer prop- 
erties of the boundary layer are concerned :° 


= R,/2 (53) 
oe Pe” 


where the formula for recovery factor is correct to 
within 0.3 per cent for 0.6 < Pr < 1.0. For this special 
case, Eqs. (19) and (20) can be integrated analytically 
(the derivatives are constant), and the results can be 
combined with Eqs. (42)—(44) to yield the following 
formulas: 


R, = V (u*/T*)w 


A*f,(0) 2 


Nur R, (54) 
(0) 
£0(0) 


Pr A** 


where A*, being constant, is determined from Eq. (31), 
which applies at the leading edge. The constant 
V 2f,(0) Bo is evaluated as 0.685 and 0.661 for 4P and 
iP, respectively, which compares well with the exact 
solution. 

In Fig. 2, A* is plotted as a function of Pr, 
together with the formula A* = Pr~'?, which would 
lead to the correct relation between the Nusselt Num- 
ber and skin-friction coefficient in Eq. (54). The re- 
covery factors of Eqs. (53) and (54) are compared in 
Fig. 7. 

The overall agreement indicated above is satisfac- 
tory. It appears, however, that the quantities in- 
volved are fairly sensitive to the choice of profile func- 
tions, which have been determined in a somewhat 
arbitrary way. 

For flows with variable surface temperature, Eq. (20) 
must be integrated numerically. This was done for the 
sample problem of Chapman and Rubesin, and the re- 
sults are compared in Fig. 8. A curve representing 
the value of the given parameter when the effect of 
surface temperature gradient is neglected has also been 
plotted. It is seen that the effect of surface temper- 
ature gradient is correctly predicted by the present 
theory, although the numerical accuracy diminishes 
with increasing Mach Number. It must be pointed 
out, however, that for this particular problem the values 
of 7,,* over the plate increase as the square of the Mach 
Number. Thus, relatively good agreement is assured 
if 7.,.* is restricted to a reasonable range. 
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SAMPLE PROBLEM 


A sample calculation has been carried out to illustrate 
the application of the present theory in the most gen- 
eral case. 
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Boundary Conditions 


The body in question was assumed to be a biconvex 
circular are airfoil of 10 per cent thickness ratio, at zero 
angle of attack. The problem is further specified by the 
following conditions: 


M. = 2.0, Pr = 0.72 
T. = —67°F., S* = 0.505 (55) 
= 2.0 — 2.0% + 0.8% 


where ¢ is the curvilinear coordinate measured from the 
jeading edge and referred to the chord length c. 

The local free-stream velocity distribution can be cal- 
culated according to Busemann’s second-order theory” 
and is written in the following form: 


fiy = + + 


to’ = + { 
where 
(y — + 2 
4(M.2 — 1)? 
(57) 


n 
=. 
= 
PS 
wad 
vr 


and 7 is the airfoil radius of curvature referred to the 
chord length. 


Method of Integration and Results 


The simple integration formula 


ly 
dx n—1 


dx),—, — (dy (x, — 2,1)? (58) 


= 
was found to be best suited to the problem at hand. The 
main part of the integration consists of the following 
operation: Given Ry? and Ro at %,, calculate the de- 
rivatives d( Ryd) ‘dé and according to Eqs. 
(19)-(24). The numerical work has been carried out 
in two different ways: (1) The entire calculation has 
been set up for the IBM Card Programmed Calculator, 
which integrates Eqs. (19) and (20) for any set of 
boundary conditions in approximately twenty minutes. 
The skin-friction and heat-transfer coefficients are ob- 
tained directly as outputs. (2) The calculation has 
been carried out by hand with the aid of a manual cal- 


culator; this requires approximately 8 hours. 
The skin-friction and heat-transfer parameters cal- 
culated for the problem specified by Eqs. (55)—(57) 


have been plotted in Figs. 9 and 10, and the important 
flow variables obtained from the solution are shown in 
Figs. 11 and 12. To establish the effect of compressi- 
bility, the low-speed and constant-density solutions are 
also indicated. Finally, the combined effect due to ex- 
ternal pressure and surface temperature gradients has 
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been indicated by plotting the local flat-plate values of 
the given parameters, as defined by Eq. (54). 


CONCLUSIONS 


The method that has been developed here is to be 
recommended both for its unusually wide range of ap- 
plication and for the ease with which numerical calcu- 
lations can be carried out. The effects of compressibil- 
ity and of pressure and surface temperature gradients 
on the characteristics of the boundary layer are ac- 
counted for simultaneously without imposing any re- 
strictions on the Prandtl Number. No special viscosity- 
temperature relation was assumed in the general de- 
velopment. In the general case, the calculation of skin- 
friction and heat-transfer parameters involves the nu- 
merical integration of two simultaneous ordinary first- 
order differential equations. Any convenient method of 
numerical integration can be utilized to reduce the cal- 
culations to a standard procedure that can be carried 
out by a computation group or by a digital calculator. 

A complete estimate of the accuracy that can be ex- 
pected from the present theory cannot be given, be- 
cause reliable solutions for comparison exist only for 
simplified flow conditions. The comparisons that can 
be made, however, indicate that the accuracy obtain- 
able is fairly good and is satisfactory for general engi- 
neering purposes. It is believed, further, that the ac- 
curacy of the method can be improved by refinement 
of the profile functions. 
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Plastic Collapse of Thin Rings’ 


CHINTSUN HWANG? 
Unwersity of Notre Dame 


SUMMARY 

When several concentrated or distributed loads are applied 
on a circular or elliptic ring and are increased proportionally, 
the ring first vields at some locations and, at the final stage, 
collapses as a mechanism. The collapse of a ring cannot take 
place until the number of yield joints developed in the ring 
exceeds its number of redundancies. In this paper we shall 
study the configuration of the collapse mechanism, the corre 
sponding collapse load of rings, and the deformations of a circular 
ring immediately before plastic collapse takes place 

For rings in which the heights of the cross section are small 
compared with their overall dimensions, the stresses are essen- 
tially due to flexure of the ring segments when the ring is under 
The stresses due to axial and 
Further- 


the action of concentrated forces. 
shear forces are relatively small and may be neglected 
more, at the instant of plastic collapse, the deformation of the 
ring as a mechanism is unlimited if we assume no strain-harden- 
ing of the material in the plastic state. In parts of the ring 
where the sections are still elastic, or partly elastic and partly 
plastic, the deformation is contained and may be neglected 
Thus the change in geometry of the ring is disregarded; the 
ring is considered as a mechanism made of perfectly rigid mem 
bers connected by vield joints that offer constant resisting 
moments. The positions of these yield joints and the corre 
sponding collapse load are found for several typical loading 


conditions. 


INTRODUCTION 


HE ELASTIC STRESS ANALYSIS of a redundant struc- 
pa under steadily increasing loads, usually a 
tedious task, is applicable only when the stresses are 
below the yield limit. The structure may be said to 
have failed when the stresses in some locations reach 
the yield stress. Another way of defining the failure 
condition depends upon the deformation of the struc- 
ture. We may say the structure fails when it yields 
completely at several locations and collapses as a 
mechanism. Thus, to study the failure of a structure 
one only needs to analyze the possible ways of total 
collapse of the structure. 

In order to produce deflections far beyond the range 
of deformations encountered in the elastic or elasto- 
plastic state in redundant structures such as portal 
frames or rings, several sections or joints must reach 
complete yielding. To satisfy the collapse condition, 
various combinations of yield joints can be assumed in 
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a structure. Each of these combinations corresponds 
to a different collapse load and to a different mode of 
failure. The actual mode of failure is the one that 
gives the smallest collapse load. 

The method of plastic collapse was first introduced 
by N. C. Kist in 1920. In 1931, Girkmann! and 
Fritsche’ used the method of plastic collapse to find 
the ultimate load of continuous beams and simple 
frames. Some later works related to this subject are 
listed in references 

A comprehensive study of the theory of plastic col- 
lapse of redundant structures was carried out by Green- 
berg and Prager.’ For a set of loads increasing pro- 
portionally, Greenberg and Prager adopted the load 
intensity at a specific location as the characteristic load 
intensity on the structure. The characteristic load 
intensity at the instant of plastic collapse was desig- 
nated as P,. For the case where the degree of redun- 
dancy of the structure was » they assigned 7 values to 
the » redundancies. A statically admissible condition 
was established in which the redundancies were stati- 
cally compatible with the external loads and no more 
than + 1 sections reached the complete yielding 
condition. The characteristic load intensity at this 
instant was designated as P,. On the ether hand, 
Greenberg and Prager also established a certain col- 
lapse mode for the structure by assuming ” + 1 yield 
hinges. The external loads that are sufficient to cause 
this particular mode of collapse were called kinemati- 
cally sufficient loads. The corresponding characteristic 
load intensity represented the magnitudes of the ex- 
ternal loads which satisfied the kinematically sufficient 
condition and was designated as P,. Greenberg and 
Prager proved that 


(1) 


From this theorem the upper and lower bounds of the 
collapse load of a redundant structure were established. 
To find the actual collapse load, it was necessary to 
adjust both statically admissible and kinematically 
sufficient conditions so that the lower and upper bounds 
would eventually coincide. In a later paper,'’ Drucker, 
Prager, and Greenberg extended the theory to more 
general types of elastoplastic materials; they also dis- 
cussed the effects of disregarding the change in geom- 
etry of the structure before and during the plastic 
collapse. 

For arches in which the axial loads may play an 
important part in causing collapse, the plastic failure 


S19 


} 


$20 JOURNAL OF THE AERONAUTICAL SCIENCES—DECEMBER, 


1953 
| 
M 10 
05 
10 20 30 40 50 60 70 80 90 
SS 
-05 Pep — 
= P. 
-1.0 
| 
p 
=p 
Fic. 1. 


The bending moment distribution along a circular ring under diametral loads. 
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of a section is defined by the axial load and the bending 
moment acting simultaneously. The collapse loads 
were analyzed in reference 11. 

In the case of thin rings, we shall study in this paper 
all possible conditions that are kinematically sufficient 
to cause plastic collapse. Then we shall try to deter- 
mine the one that corresponds to the minimum char- 
acteristic load intensity. This minimum value repre- 
sents the actual plastic collapse load intensity of the 
structure. 


PLASTIC COLLAPSE OF A THIN CIRCULAR RING UNDER 
THE ACTION OF TWo-DIAMETRAL FORCES 


We start with the simple case of a thin circular ring 
under the action of two-diametral forces. It is as- 
sumed that the cross section of the ring has two axes 
of symmetry perpendicular to each other and that one 
of the axes of symmetry is parallel with the ring axis. 
The height of the cross section is assumed to be small 
compared with the overall dimensions of the ring. 
When external loads are acting on the ring, the stresses 
in the ring sections are essentially due to flexure of the 
ring segments. The stresses due to axial and shear 
forces are relatively small and may be neglected. The 
material of the ring is assumed to be perfectly elasto- 
plastic, with identical stress-strain diagrams for both 
tension and compression. 

The two equal and opposite forces are applied at the 
points A and C as shown in Fig. 1. If the bending 
moment at A is designated by 1/4, the bending mo- 
ment M at any point (R, 8) is 

M = M, — (1/2)PR sin 6 (O<@0< 7) (2) 
By means of the energy method we find 


Ma = (3) 


when the ring is perfectly elastic. 
Let us increase gradually the intensity of the loads P. 
For 


P = P, = xo,I/(Re) 


_the yield stress o, is reached in the outer fibers of 


The corresponding bending mo- 
If we 


sections A and C. 
ment is ¢,J/c, where //c is the section modulus. 
increase the intensity of the loads P beyond the value 
P;, then for a value P = P» the entire sections A and 
C become plastic. When parts of the ring are plastic, 
Eq. (3) is no longer valid because it is derived from the 
elastic energy theory. Yet for P = P2, the yield region 
is small compared to the part of the ring which is still 
elastic, and therefore the ring as a whole can be con- 
sidered as elastic with a plastic hinge at each of the 
cross sections A and C. Hence, from Eq. (3), we 
find 

P, = (4) 


where S is the statical moment of the area of the upper 
or lower half of the cross section with respect to the 
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centroidal axis. Even after the value P: has been 
reached, the ring can take more loads, because up to 
now only two of its sections have yielded. 

If plastic hinges also appear at the cross sections at 
B and D, then the collapse load P, is reached. The 
value of P, may be found by setting JJ = \/, at both 
6 = Oand = The result is 


P, = 8So,/R = 4M,/R (5) 


The moment distributions along the are AB of a ring 
with a rectangular cross section under the influence of 
P,, P2, and P, are shown in Fig. 1. At sections A, B, 
C, and D, Fig. 1, four yield hinges are formed at the 
moment of plastic collapse. 

The collapse load P, may also be found by consider- 
ing the energy dissipated in the ring. The work done 
by the equal and opposite forces P, in moving a dis- 
tance d (Fig. 2) is —2P,.RAa, while the amount of the 
energy dissipated in the four yield hinges at A, B, C, 
and D for a change in angle Aa is —4.17,(2Aa). Equat- 
ing the work done by the forces P, to the dissipated 
energy, we obtain 


P, = 4M,/R 


The simplicity in the solution of the collapse load is 
due to the fact that the positions of the four yield 
hinges can be assumed in view of the symmetry of 
both the structure and the external loads. When the 
position of the yield hinges is not obvious, we resort to 
a more refined method in order to find the collapse loads, 
and this computation will be shown later. 

When P < P2, the deflection 6 = 2d, Fig. 2, along 
the x-axis is 


6 = 0.1488PR*/(ET) (6) 


In the case of a perfectly elastic structure, the de- 
flection is found by the well-known Castigliano’s 
theorem, 


5 = 0U/aP 


When a structure reaches the elastoplastic state, the 
flexural strain energy can be represented as 


a M 
U= If Mdads = [Ma — a dM) ds 
s 0 s 0 


where a = a(J/) is the angle of bending per unit length 
and s is the length measured along the bending segment; 


the term 
M 
adM 
0 


is called the complementary energy. We differentiate 
the strain energy UL’ with respect to an external load P 


and obtain 
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= 0U/OM-OM/OP 


the section, the strain energy in the ring is, approxi- 


= Qa d as n/2 9 
~ a as = 


Applying Eq. (7), we obtain 
where 6 is the deflection in the direction along which 


the load P is applied. In working out Eq. (7), the J PRsin a) Resin 0 (9) 
only approximation involved is EI Jo 2 


M(Oa/OM) ~ a When P > P,, Eq. (9) gives the ring deformation im- 


mediately before the complete plastic collapse takes 
place, 


which is reasonable if the structure is within the elasto- 
plastic state and if the moment-angular deflection curve 
has a sharp defined turn when the yield moment /, is 
reached. 

For the circular ring under the diametral load P, 
(P, < P < P,), because the ring segments are either 
in the elastic state or contained by the elastic core of 


6. = 1.1416./,R2/(ED 


Fig. 2 is a plot of the deformation 6 of the ring of 
rectangular cross section versus increasing diametral 
loads P. 


When the loads P are applied, the circular ring deforms into an oval shape. We consider the displacement eé 
of the section B when the diametral forces are applied along the x-axis (Fig. 3). Assume a pair of virtual loads 
Q applied along the y-axis toward the origin at the sections B and D. Then, for P < P2, we have 


M = ((P/m) — (1/2)P sin — (Q/mr) + (1/2)0 cos [(0<6< (1/2)r] (10) 


Applying Eq. (8) to the ring under the additional virtual loads Q, we find the deflection e at section B to be 


us 


] ) ] 1 
e= {12° (EI) | ( P sin 6 — ¢ + = Q cos cos — = 0.06831PR*/(EI) (11) 


In the loading range P, < P < P,, we have 
‘M = M, — (1/2)PRsin — (1/2)QR(1 — cos @) (0 < 6 < (1/2)r] (12) 


The corresponding deflection e of the section B along the y-axis is 


ah, PR sin 6 —_ 9 GRU cos) 0) R? = 


( —0.5708M, PR)R: (ED) (13) 


The ratios 6/(2R) and eR versus P/(.\/,/R) are computed for various values of R/k, k being the length of the 
sides of the square cross section. The result is plotted in Fig. 3. In these calculations ¢,// is assumed to be 
1/1,000. It is interesting to note that, even for R/k = 20, the maximum change in length of the diameter in the 
y-direction is 2.58 per cent, while the maximum change in length of the diameter in the x-direction is 3.42 per 
cent. In view of the fact that the deformation is small before plastic collapse takes place, the bending moment 
along the ring is not affected to any appreciable degree by this deformation. 


PLASTIC COLLAPSE OF THIN ELLIPTIC RINGS 
UNDER CONCENTRATED FORCES 


Burke in the form of working charts.'‘ The stress 
concentration at the points of application of the con- 
centrated forces is neglected, and the possibility of 


Imagine a thin elliptic ring acted upon by concen- 
trated forces that are symmetrical with respect to the 
major axis of the ring, as shown in Fig. 4. In the 
elastic state the redundant bending moment and axial 
force may be determined from energy considerations. 
A complete set of such redundant data was given by 


elastic buckling of the ring will not be studied here. 
When the material is stressed beyond its elastic limit 
and when the external forces are steadily increased, 
some sections in the elliptic ring yield completely. 
These sections behave as yield joints with a constant 
resisting bending moment ./,. Four such yield joints 


pares tie 
822 | 
mately, 
= 
are } 
four: 
Beca 
colla 
to tk 
othe 
by t 
in al 
T 
com 
sect 
(Ba 
enin 
mon 
and 
3 In t 
are 
the 
the 
A 
ring 
| whe 
two 
| 
= | 


TOXi- 


(8) 


x of 
‘tral 


nt e 
ads 


10) 


PLASTIC COLLAPSE OF THIN RINGS 823 
(c) (c) (¢) (c ) 
R & 
(2) 42) 
30 | / / J) | 
| / ; 
Q 
(1) (1) 
: : a 
| 
| Q 
| ff | 
0 004 008 ol2 O16 020 024 028 
R 
Fic. 3. The diametral deformation of a circular ring along the y-axis. 


are necessary to cause the ring to degenerate into a 
four-bar mechanism at the instant of plastic collapse. 
Because of the symmetry of the forces with respect to 
the major axis of the elliptic ring, we may expect the 
collapsing mechanism to be symmetrical with respect 
to the major axis. Hence, two of the yield joints should 
be located at the terminals of the major axis. The 
other two yield hinges, which are also located sym- 
metrically with respect to the major axis, must be found 
by the condition that the collapse load is the smallest 
in all the kinematically sufficient loads. 

To take care of the local stress concentration, it is a 
common engineering practice to strengthen the cross 
sections in the neighborhood of the concentrated forces 
(B and F in Fig. 4). In the elastic state, this strength- 
ening will influence the magnitude of the redundant 
moment, as well as that of the redundant axial force, 
and change the moment distribution along the ring. 
In the case of plastic collapse, as long as the sections 
are not changed where complete yielding takes place, 
the strengthening of parts of the ring will not influence 
the collapse load and the positions of the yield joints. 

Adopting the minor and major axes of the elliptic 
ting as the Y and J axes, the centroidal line of the 
ting may be defined by 


X = bdbecosy, 


asin 


where y is a parameter varying from —7 to mr. The 
two points of application B and F of the concentrated 
forces are designated by 


Xp = —Np cos do 


= Yr = asin 


where @» is the eccentric angle corresponding to points 
Band F. Under the action of the concentrated forces 
and when plastic collapse is imminent, two of the four 
yield hinges are formed at the terminals A and D of the 
major axis (Fig. 4). The other two yield sections C 
and £, whose positions are still unknown, may be 
designated by 


= 
Vr =asing 


Yo 


where ¢ is a quantity to be determined. 
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The collapse mechanism of an elliptic ring. 
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venience in calculation we shall adopt a new coordinate 
system (x,y) in which the y-axis coincides with the 
}-axis while the x-axis passes through points C and 
(Fig. 4). 

When the four sections A, C, D, and E yield com- 
pletely, the ring becomes a mechanism. The coordi- 
nates of the four yield sections (x4, va, .. .) and those of 
the points of application of the concentrated forces 
(xp, Va, ...) are indicated in Fig. 4. 

From geometrical considerations we have 


Xc 


Vb? cos? @+ + sin 


cos a = (14) 


Vb? cos? @ + + sin 


(15) 


sin a = 


At the instant of plastic collapse, the deformed ring 
remains symmetrical with respect to the y-axis. Thus 
a small change Aa in the angle a will be sufficient to 
define the corresponding change in positions of the 
hinges, as well as those of the points of application of 
the concentrated forces. Thus we have 


Axe = —a(1 + sin ¢)Aa (16) 
Ayp = —bcos ¢ Aa (17) 


From geometrical considerations we also have 
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For an increment of the angle Aa, the deformation of 
the point B is 


(1 + sin ¢)(1 — sin ¢o) 
a Aa 


Ax = — 9 
~ 1 — sing (20) 
along the x-direction and 
(cos @ — cos go)(1 + sin Aw (21) 


1 — sing 


along the y-direction. 

Corresponding to a change of angle Aa of the mech- 
anism at the instant of plastic collapse, work is done 
by the external forces because their points of applica- 
tion are moved. This work is transformed into dissi- 
pated energy in the four yield hinges. From this 
postulation we can find the magnitudes of the collapse 
loads for various positions of the yield hinges C and E, 
The true position of the yield hinges will be the one that 
corresponds to the minimum collapse loads. For a 
change in the angle Aa, the energy dissipated in the 
four hinges of the ring is 


A(D.E.) = 4M,(Aa + 
= 8M,Aa/(1 — sin ¢) (22 


The corresponding increment in the work done by the 
forces P and Q is 
‘AW = 2P(—Ayp) + 2P(Ayz) — 


1+ sing Substituting values for Ayp, Ayg, and Axg from Eqs. 
aie ( — sing ae (18) (17), (21), and (20) and using the notation 
i = Qa/(Pb) 
1 + sin *) m = Qa/( 
Aya = A 19 
” ald (; — sing a (19) the work done by the forces P and Q is reduced to 
AW = 2PbAa — cos go — Cos sin + m(1 + sin (23) 
1 — sing 


Equating the work done by the forces to the corresponding dissipated energy, we obtain 


4M 1 
( b 2 cos ¢ — cos do — cos gy sin d + m(1 + sin d)(1 — sin gd») 


(24) 


In Eq. (24), m and @q» are given quantities, while ¢ is an unknown quantity to be determined from the condition 


that the collapse loads P should be a minimum. 
Applying this condition to Eq. (24), we find 


—2 sin — cos cos + m(1 — sin go) cos 


For minimum collapse loads P, we require that 0P/0¢ = 0. 


2 cos d — cos do — cos go sin d + m(1 + sin d)(1 — sin do) 


which can be reduced to 


tan @ = (1/2)[m(1 — sin go) — cos go] (25) 


Eq. (25) determines the position of yield hinges C and E (Fig. 4) for given values of m and @. Thus, for any 
thin elliptic ring acted upon by a set of loads as shown in Fig. 4 and as long as the ratio m = Qa/(Pb) and the 
position of the points of application of the concentrated forces are given, we can find the position of the corre- 
sponding yield hinges. Eq. (25) is also plotted in Fig. 5 for various values of m. 
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Once the value ¢ is determined, we may calculate the collapse load P, from Eq. (24). 


A chart for the determination of the positions of yield sections of an elliptic ring under symmetrical concentrated forces. 


By substituting the 


value of ¢ from Eq. (25) into Eq. (24) and after some simplification, we find the collapse load to be 


1 
(26) 


P. = (2) 
b m(1 — sin @o) — cos do + V4 + [m(1 — sin — cos 


To illustrate the results of our analysis, we shall calculate the collapse loads, the position of the plastic hinges, 


and the moment distribution for a thin elliptic ring. 


EXAMPLE 


Let us assume that the ratio between the semimajor 
and semiminor axes to be a/b = 3/2 and that the forces 
P and Q are acting at sections B and F, which are 
located by radius vectors making an angle of 45° with 


the minor axis. Thus we have 


Vp/Xp = 1 = asin cos dp) 
which gives a value of @ = 33.7°. If Q = (1/2)P, we 
have 


m = Qa/(Pb) = 3/4 


By aid of reference 14, the bending moment diagram 
of the elliptic ring when one of the sections (section D) 
yields completely is plotted in Fig. 6. 
sponding value of P is found to be 


The corre- 


1.51.,/b 


When complete plastic collapse takes place, the 
positions of the yield hinges are found by substituting 
the values of m and @» in Eq. (25), or by using data in 
Fig. 5, the result is ¢ = —14.0°. Then, by making 
use of Eq. (26), we find 

P, = 2.56M,/b 


The bending moment diagram for a half span of the 
ring corresponding to P = P, is also plotted in Fig. 6. 


CONCLUSION 


In this paper only the case of plastic collapse of thin 
elliptic rings under concentrated forces that are sym- 
metrical to the major axis was studied. Actually there 
is no limitation to the shape of the ring and the arrange- 


| 
(21) 
lech- 
done 
lica- 
lissi- 
this 
apse 
d F, 
that 
ra | 
the 
| 
-454 
= 
4) 
on 
0. | 


826 JOURNAL OF THE AERONAUTI 


CAL 


SCIENCES—DECEMBER, 1953 


ate 
<— 


QO -45 -30 O 


90 


I5 30 45 //60 75 
Wa 


\ 


Fic. 6. 


ment of the concentrated forces. For a thin oval- 
shaped ring under several concentrated forces, its 
plastic collapse condition may be studied in a manner 
similar to the analysis of the elliptic rings. To analyze 
the oval-shaped ring, we shall need again an analytical 
expression for its centroidal axis, or the process may 
be performed graphically. The positions of the yield 
hinges and the collapse loads are computed by mini- 
mizing the loads corresponding to all kinematically 
sufficient conditions. However, the oval-shaped ring 
analysis will be complicated if the collapse configura- 
tion has no axis of symmetry. 

Several basic assumptions have been made in this 
paper, and these are now reviewed. Since the material 
was assumed to be perfectly elastoplastic, the effect of 
strain-hardening was ignored. Also neglected was the 
effect of the expansion of the plastic yield region in the 
neighborhood of the yield sections. Symonds and 
Neal discussed these two factors in reference 8, and 
they believed that the errors introduced by these two 
assumptions tend to cancel each other. It is thus 
expected that the collapse loads derived under these 
assumptions are fairly correct. 

Another item that should be considered is the 
deformation of the ring before complete collapse takes 
place. If it is excessive, the distribution of the bending 
moment may be influenced substantially, and an 
additional error is involved in the collapse loads. In 
the case of a circular ring we investigated the deforma- 
tion immediately before plastic collapse takes place 


The bending moment distribution along an elliptic ring under concentrated forces. 


and found the deformation to be small; hence, the 
error involved in determining the collapse load is 
negligible. 
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Effects of Wing Twist on the Response of an 
Airplane Encountering a Sharp-Edged Gust’ 


GORDON C. K. YEH? ann JOHANN MARTINEK? 


Reed Research Inc. 


SUMMARY 


An analysis is made of the gust response of a simplified model 
of an airplane having the degrees of freedom of vertical motion 
and wing twist. A two-dimensional, incompressible unsteady 
flow isassumed. Under the further assumption of linear relation 
between the wing twist and the vertical acceleration, the equation 
of motion is set up as a single integrodifferential equation and is 
solved by the method of Laplace transformation with the use of 
the Faltung theorem. The acceleration responses and ailevia- 
tion factors are then computed and plotted for a systematic vari- 
ation of two dimensionless parameters \ and cR, representing the 
mass ratio and the wing torsional flexibility, respectively. The 
resulting curves for conventional swept-back wing aircraft show 
that: 

(a) The wing elasticity produces considerable alleviation effects, 
especially for small values of X. 

(b) It has negligible influence on the time at which the maxi- 
mum acceleration is attained. 

(c) The rate of decrease of the alleviation factor with increase 
in cR diminishes as cR increases and tends to be constant for large 
values of X. 


SYMBOLS 

W = aircraft weight 

M = aircraft mass 

b = wing span 

C = wing mean aerodynamic chord 

S = wing area = be 

L = aerodynamic lift per unit span, positive upward 

L, = lift per unit span due to gust only 

| = lift per unit span due to the motion of the airfoil 
only 

) ae = L,, due to the apparent mass only 

Lay = L,, due to the quasi-steady and wake effects 
only 

p = air density 

U = aircraft forward velocity 

t = time variable 

s = dimensionless variable = (2/c)Ut 

o = variable of integration in place of s in Duhamel 
integrals 

£ = Laplacian operator 

p = variable in Laplace transforms 

z = vertical displacement of aircraft c.g., positive 
downward 

6 = angular displacement of the wing from its ini- 
tial position, positive clockwise 

a = angle of attack of the wing 

ay = initial (steady) angle of attack of the wing be- 


fore entry into the gust 
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m = slope of the lift coefficient curve for the aircraft 
= dC, /da 

= dimensionless parameter = 1//p(c/2)Sm 

R = ratio between @(s) and s"(s), as defined in Eq. 
(1) 

los, dso, bas = chordwise distances shown in Fig. | in terms 
of the half-chord 

w = vertical gust velocity 

W(s) = K4arman-Sears function for the lift following 


penetration of a sharp-edged gust in incom- 
pressible flow 


(5) = Wagner function for the lift following a sudden 
change in angle of attack in incompressible 
flow 


43, Q2, 4), 2g = coefficients of a cubic equation Eq. (A6) 


by, be, bs = roots of the cubic equation Eq. (A6) 
N(p) = a function of p, as defined in Eq. (A2) 
yp) = a function of p, as defined in Eq. (A3) 
ad p) = a function of p, as defined in Eq. (A5) 
F(s) = acceleration response 


K alleviation factor 
oY = wing twist alleviation factor, as defined in 
Eq. (36) 
Dots denote derivatives with respect to ¢; prime marks denote 
derivatives with respect to s. Other symbols will be defined as 


they are introduced. 


INTRODUCTION 


WwW" THE INCREASE IN SIZE AND SPEED of present- 
day aircraft, the increase in the response of an 
airplane to a gust has caused critical conditions in all 
phases of aircraft design. Therefore, the number of 
papers dealing with the prediction of wing loads on an 
airplane due to a vertical gust has grown significantly. 
Many papers have treated the airplane as a rigid body 
and, in so doing, have dealt with either the degree of 
freedom of vertical motion alone or with the degrees 
of freedom of vertical motion and pitch, such as the 
papers mentioned in references 1 and 2. Of greater 
concern in the consideration of gust penetration is the 
influence of wing flexibility on airplane response, par- 
ticularly in the case of a swept-back wing, for which 
twist effects and the change of angle of attack due to 
bending may result in considerable alleviation of 
gust load. This has given rise to many recent papers 
that treat the airplane as an elastic body. (See, for ex- 
ample, references 1, 2, and 3.) Invariably, these 
elastic-body analyses have the disadvantage of being 
laborious in computation and are not suitable for making 
trend studies. 
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Fic. 1, Configuration and symbols for the model wing. 


In the present report, the case of the gust penetration 
of an airplane having the degrees of freedom of vertical 
motion and wing twist is considered. While several 
papers have dealt with the degrees of freedom of verti- 
cal motion and wing bending (see, for example, refer- 
ence 4) or the degrees of freedom of vertical motion 
and wing bending and twist (see, for example, reference 
5), no work has yet presented a systematic set of results 
of acceleration responses for the degrees of freedom of 
vertical motion and wing twist alone. This paper has 
the objectives of: (a) establishing some of the basic 
parameters that are involved when wing twist flexi- 
bility is included; (b) developing a method of solution 
which is fairly well suited for trend studies without ex- 
cessive computation time; (c) presenting curves show- 
ing the basic parameter variations that will provide in- 
sight into the nature of the effects on the acceleration 
response of an airplane traversing a gust due to the 
wing twist alone. 


THE SIMPLIFIED MODEL AND ASSUMPTIONS 


In order to present our method of analysis for gust 
loads on airplane wings without undue complication, 
the airplane will first be replaced by a simplified model. 
The analysis of this simplified model retains the essential 
features of the analysis of the actual airplane with wing- 
twist flexibility and is much easier to follow. Other as- 
sumptions used in the analysis are also stated here. 

The model consists of a rigid and massless rectangular 
wing (Fig. 1), which is connected by a rigid post 
(fixed to the fuselage and hinged to the wing) and a 
torsion spring to a fuselage of finite mass and infinite 
pitching moment of inertia. The point of connection 
of the rigid post and the torsion spring to the wing cor- 
responds to the elastic center of an actual wing. The 
flow about the wing is assumed to be two dimensional 
and incompressible. It is assumed that there are no 
aerodynamic forces acting on the fuselage. 

The rigid post signifies that the bending stiffness of 
the wing is so high that it may be assumed to be in- 
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finity, while the torsion spring permits the wing to ro- 
tate in the pitching direction, and simulates the tor- 
sional stiffness of an actual wing. The present model 
assumes that this torsional stiffness is such that the 
wing twist @(s) is directly proportional to the vertical 
acceleration of the airplane 2”(s) as a function of the 
nondimensional parameter s: 


= Rs"(s) (1) 


This follows froin the assumption that the wing is mass- 
less, which implies that the moment of inertia of the 
wing about the elastic axis is zero. This assumption 
does not apply to transient stress and flutter analyses 
but is considered to be reasonable for acceleration re- 
sponse analysis. 

The model is assumed to be in steady horizontal flight 
at velocity U before entry into a vertical sharp-edged 
gust of velocity w and uniform spanwise distribution. 
Even after entry into the gust, the horizontal compo- 
nent of velocity is assumed to remain constant at the 
value U. The model is thus seen to possess two degrees 
of freedom—a rigid body degree of freedom for transla- 
tion in the vertical direction and a vibratory degree of 
freedom for elastic angular displacement of the wing 
relative to the fuselage. Through the assumption Eq. 
(1), however, the problem is reduced to that of one de- 
gree of freedom. Only motion in the plane of sym- 
metry is possible since the gust is symmetrical. 

The dimensional properties of the coordinate system 
of the wing are shown in Fig. 1. ag is the steady angle 
of attack of the wing before entry into the gust. 2 and 
6 are coordinates that measure the vertical and pitch- 
ing displacements of the wing from its initial position, 

‘with respect to fixed axes, and are positive in the direc- 
tions shown. The quantities 1/5, Js, etc., define the 
chordwise distances shown in terms of the half-chord. 

The angles ap and @ are assumed to be small, and the 
velocities w and % are assumed to be small compared 
with U. On the basis of these assumptions, the follow- 
ing approximations may be considered valid: 

(1) The magnitude of the relative wind velocity is at 
all times equal to U. 

(2) The component of the aerodynamic force in the 
z-direction is equal to the lift. 

(3) The aerodynamic moment about the elastic center 
is equal to the product of the lift and the distance from 
the chordwise location of the lift vector to the elastic 
center. 

A further assumption with regard to the aerody- 
namic forces acting on the wing is that the part of appar- 
ent mass forces due to the wing twist and the part of 
quasi-steady and wake forces due to /7; ~ 0 are so small 
that they may be neglected when compared with other 
aerodynamic forces. This assumption reduces the 
amount of numerical computations considerably, and 
yet Appendix (B) shows that they are well justified 
within the range of conventional airplanes. 
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EXPRESSIONS FOR THE AERODYNAMIC FORCES 


In the present problem the forcing function on the 
system is derived entirely from the unsteady aerody- 
namic forces acting on the wing, since the gravity force 
and the aerodynamic forces acting in the initial steady 
flight condition are in equilibrium and therefore cancel 
out. Because of the simplifying assumption, Eq. (1), it is 
sufficient for us to consider only the disturbance lifts 
without considering the disturbance moments. In the 
determination of the general expressions for these dis- 
turbance lifts, we make use of an idealized theory for 
two-dimensional flow about a thin airfoil in unsteady 
motion in a nonviscous, incompressible fluid. With 
the additional assumption that the angle of attack, due 
to any cause, is at all times small, the theory leads to 
aerodynamic force expressions that are linear functions 
of the angle of attack, the velocities, and the acceler- 
ations of the airfoil. Thus the principle of superposi- 
tion applies in the present case. 

We express the net disturbance lift per unit span L 
acting on the airfoil as follows: 


(2) 


where L, and L,, are, respectively, the lift per unit 
span due to the gust only and the lift per unit span due 
to the resulting disturbed motion of the airfoil only. 
L,, can be further decomposed as follows: 


Limp + | (3) 


where /.,,, is the lift on the airfoil due to the quasi- 
steady and wake effects and L,,, is the apparent mass 
lift. 

In reference 1, it is shown that, for a model wing 
with an arbitrary location of the elastic center, as shown 
in Fig. 1, these component lifts take the following ex- 
pressions: 


1 *s 
L,(s) = w(a)V'(s — «) da (4) 
0 


Lnp(s) = (s — a) X 
Jo 
4- + | da (5) 


Lme(S) = (1/4) + 
6’(s) + (6) 


where V(s) is the von Karman-Sears function for the 
lift following penetration of a sharp-edged gust in in- 
compressible flow® and #(s) is the Wagner function for 
the lift following a sudden change in angle of attack in 
incompressible flow.* 7 The exact functions ¥(s) and 
#(s) cannot be expressed in simple analytical form. 
Satisfactory approximations, however, are given by the 
following expressions :* 


V(s) = 1 — — 0.500e-* (7) 
= 1 — — (8) 


THE EQUATION OF MOTION AND ITs SOLUTION 


Let us consider the disturbed motion of the simplified 
model. Since the pitching moment of inertia of the 
model is infinite, only its vertical motion needs to be 
considered. This disturbed vertical motion is governed 
by the following differential equation: 


M3(t) = —bL(H (9) 
By using the parameter 
s = (2/c)Ut (10) 


which is the distance traveled in half-chords, Eq. (9) 
becomes 


MU*(2/c)*2"(s) 


—bL(s) 
—bL,(s) — bLm,(s) — bLma(s) (11) 


For a sharp-edged gust 
w(s) = w = constant (12) 


L,(s) can be expressed from Eq. (4) as follows: 
1 s 

pUeme — «)de 
0 


= (1/2)pUcmw[¥(s) — ¥(0)] (13) 
Although, in general, 
ls # 0 (14) 
the forces resulting from this are small so that the term 
1,0” (s) in Eq. (5) can be neglected. Therefore, we have 
the following expression for L»,(s): 
1 s 
Lny(s) = @(s — a) X 
0 
9 
ie | da (15) 
From the assumption that the part of apparent mass 
forces due to the wing twist is so small that it may be 


neglected compared with other aerodynamic forces, 
Lima(S) can be derived from Eq. (6) as follows: 


Ling(S) = (1/2) pU*cm[(1/c)2"(s)] (16) 


Substituting Eqs. (13), (15), and (16) into Eq. (11), 
we have 


MU? 2"(s) = eUSmr{¥(s) — ¥(0)] — 


s 9 
peursm — a) + ds — 
0 
5 pl ‘sm| 9 | (17) 
where the relation 
S = bc (18) 


has been used. 
Now with the assumption Eq. (1), we reduce Eq. 
(17) to the form 
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MU? (2)' =" = pU2Sm [¥(s) — ¥(0)] — 
0 
2 c 


which, with the notation 
- 2W/S M (20) 
pgcm p(c/2) Sm 
becomes 
= — ¥(0)] +f — a) X 
0 
2 1 
R2/""(o) + (9) do + (21) 


After rearranging terms and noting that, from Eq. 


(7), 
v(0) = 0 (22) 


Eq. (21) takes the final form 


cw 
= (4d + 1)3"(s) = FP + 


0 


cR f — a) da (23) 
0 


This is the integrodifferential equation of motion for the 
simplified model under the assumed conditions. 

It is observed that the form of Eq. (23) is convenient 
for solution by the method of Laplace transformation 
wherein, according to the Faltung theorem, the trans- 
form of a Duhamel integral is the product of the trans- 
forms of the two functions involved. (See, for example, 
reference 10, page 37.) 

Using the following notations and rules (reference 
10, page 294): 


~ 
° 


ACCELERATION RAT 10 Fis) 


A= 13.33 


° 


60 n= 5.34 


oOo 


10 20 30 40 50 


Fic. 2. Nondimensional acceleration responses for cR = 0. 
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(1) £{8(s)} = 
(2) = | 
(3) £{2"s)} = | 
(4) = = — 0) | 


(5) £{AG(s) + BH(s)} = AL{G(s)} + 
Bg{H(s)} = AG(p) + 


(where A and B are constants) 


(6) ef F\(s — Fo(o) do! = 
£{F,(s)}-2{ F(s)} = F\(p)*F(p) 


we obtain the Laplace transform of Eq. (23) as follows: 


(40 + 1)8"(p) = (cw/U)W(p) + 


[(2 + cRp)2"(p) — cR2"(+ 0))®(p) (25) 


Before entry into the gust the model is assumed to 
be in steady horizontal flight in which the gravity force 
and the aerodynamic forces are in equilibrium. There- 
fore, the following initial conditions hold: 


2(+0) = s’(+0) = 2"(+0) =0 (26) 


Using the initial condition 3”(+0) = 0 and rearranging 
the terms, we obtain an expression for 2”(p) from Eq. 
(25) as follows: 


) 
(4. + 1) + (2 + cRp)B(p) 


From the definitions of 2”(s) and 2(t), we have the 
relation 


a) = 02) = 
4Uw \ 


—3"(p) = (27) 


+ 1) + (2 + cRp)B(p): 
(28) 


Also, from Eq. (2.11) of reference 2, the maximum ac- 
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VS: 


SHARP-EDGED 


celeration given by the sharp-edged gust formula (for a 
rigid airplane in a steady flow) for a gust with the same 
velocity intensity w as what the model encounters is 


= —(1/A) (Uw/c) (29) 


max 
Hence, from the definition of acceleration response we 
have 


VA 


F(s) = acceleration response = >—. = 


GUST RESPONSE 


ation factor (see, for example, reference 2, page 10): 
K = (31) 


To evaluate F(s), we first proceed to obtain the ex- 
pressions for ¥(p) and ®(p) by applying the rules of 
Laplace transformation” to Eqs. (7) and (8) as fol- 
lows: 


= 0.565p + 0.130 
4dW(p) (30) p(p + 0.130) (p + 1) | 
(32 
+ 1) + + Bip) — + 0.2808p + 
whose maximum value can be regarded as an allevi- " p® + 0.3455p? + 0.01365p 
Substituting Eqs. (32) into Eq. (30), we have 
Fis) JAN(0.565p + 0.130) (p? + 0.3455p + 0.01365) | (33) 
Ss = oh 
sp? + + arp + ao) (p + 0.130) (p + 
where 
as = 1 + 44 + 0.500cR ] 
a2 = 1.3455 + 1.3820A + 0.2808cR (34) 
a, = 0.57525 + 0.0546 + 0.01365cR ’ 
a) = 0.0273 


For cases where a; is not zero (which include all cases where cR > 0) Eq. (33) may then be written in final form as 


2.260 


a3 


F(s) = 


(p + 0.130)(p + 1)( 
3 a3 a3 


(p + 0.2301) (p? + 0.3455p + 0.01365) _ j 


(35) 


In Appendix (A) it is shown that, for any given combination of \ and cR, the inverse transform in Eq. (35) can 


be easily expressed as a function of s according to the Heaviside expansion theorem. 


Hence, a response curve 


F(s) can be plotted for each combination of \ and cR.. The maximum value of F(s) is the alleviation factor K for 


the specific combination of \ and cR considered. 


RESULTS AND CONCLUSIONS 


Numerical computations have been carried out for 
20 combinations of \ and cR resulting from five values 
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of \ and four values of cR representing conventional 
The resulting acceleration 
From the 


swept-back wing aircraft. 


response curves are shown in Figs. 2-5. 
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Fic. 7. Curves of alleviation factor A vs. cR. 


maximum values of these curves, the curves of allevi- 
ation factor K are plotted in Figs. 6 and 7 against A and 
cR, respectively. Finally, based upon the data in 
Fig. 6, the curves of a so-called “wing twist alleviation 
factor” y are plotted in Fig.8. This factor is defined as: 


= (36) 


and is a measure of alleviation effect due to the wing 
twist only. For an airplane without wing twist cR = 
Oand y = 1. 

It will be noticed that for the case of cR = 0, which 
represents a rigid airplane without the wing-twist de- 
gree of freedom, the curves in Fig. 2 check eminently 
with those presented in Fig. 2.7 of reference 2. 

In conclusion, the curves of acceleration responses and 
alleviation factors for a simplified model of an airplane 
encountering a sharp-edged gust have been obtained 
for a systematic variation of the dimensionless param- 
eters \ and cR. represents the configuration and the 
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Fic. 8. Curves of wing twist alleviation factor. 


mass of the airplane, and air density and cR represents 
the torsional flexibility of the wing. [See Appendix 
(C) for the general meaning and the estimate of cR.] 
The resulting curves show that the wing flexibility of 
swept-back wings can produce considerable alleviation 
effects, especially for small values of \. Figs. 2 to 5 
show that wing elasticity tends to flatten out the peaks 
of the acceleration responses, although it does not seem 
to have much influence on the time at which the maxi- 
mum acceleration is reached. Fig. 7 shows that the 
rate of ‘decrease of A with increase in cR is greater for 
small values of cR than for large values of cR. This 
is particularly evident for small values of \. For 
large values of \, this rate of change tends to be con- 
stant. 

For a certain range of Mach Number, .\/,, where the 
Prandtl-Glauert correction may be assumed _ valid 
these resulting curves could be considered to apply 
equally to the compressible flow case with the parameter 
\ redefined as follows: 


d= MV1 — M?/p(c/2)Sm (37) 


For higher Mach Numbers and for transonic flow, more 
rigorous analyses must be used. 

It is obvious that the analysis can be easily extended 
to study the responses due to gusts of shapes other 
than sharp-edged and also to the responses due to gusts 
in supersonic flow if the corresponding functions of 
w(s), W(s), and #(s) are used. 

Further research on the influence of pitching degree 
of freedom and modification of lift functions for various 
wing plan forms, as well as compressibility effects, 
would be worth while in order to identify the most im- 
portant parameters in realistic cases. 


APPENDIX (A).—EVALUATION OF F(s) FROM Eq. (35) 


Let us express Eq. (35) as follows: 


as 


2.260A 


F(s) = 


832 


} 


where 


q(P) 


Then 
roots 


wher 


Eq 


and | 


Fc 


Eq. | 
N(- 


Si 
(A6) 


who: 


acce 


wak 


7 
= 
= 
= 
AES, 
q 
Ir 
th 
= 


SHARP-EDGED GUST RESPONSE 833 


where 


N(p) = (p + 0.2301) (p? + 0.3455p + 0.01365) = p* + 0.5756p? + 0.09315p + 0.003141 = (A2) 


g(p) = G + + + (p + 0.130) (p + 1) = (p — by) — be) (p — 43) (p + 0.130) (p + 1) 
(A3) 


Then, according to the Heaviside expansion theorem (see, for example, reference 10, page 44), if q(p) has distinct 
roots (real or complex), Eq. (A1) can be expressed in the form 


a3 N(b2) N(bs) , N(—0.130) M(—1) _, 
where 
= — bn) (A5) 
Eq. (A3) implies that },, b:, and b; are three roots of the following cubic equation: 
+ (a2/as)p? + (a:/a3)p + (ao/as) = 0 (A6) 
and Eq. (A5) means that the following expressions hold: 
= = —(b, — be) (bs — bi) + 0.130) + 1) 
= = —(b, — bs) (bz — bs) (b2 + 0.130) + 1) 
p — be p = bs 
_ gr) 
g3(b3) = i b = —(b; — b;) (be — b3) (bs + 0.130) (b3 + 1) (A7) 
2? — = bs 
gs(—0.130) = E 0.870(b; + 0.130) (b2 + 0.130) (bs + 0.130) 
qs8(—1) = Fea = +0.870(b: + 1) (b2 + 1) (bs + 1) 
For a complex value of /, say 
p= —b + ai (A8) 


Eq. (A2) becomes 
N(—5 + ai) = —[b(0.09315 — 2b?) — (6? — a*) (0.5756 — 36) — 0.003141] + 
(3b? — a? — 1.15126 + 0.09315)4i (A9) 


Since complex roots always appear in pairs, there will be at least one real root, say b;, for the cubic equation Eq. 
(A6). When d, is known, the equation can be immediately reduced to a quadratic equation, 


1 a a a 
pt + ap + — 4 4%) (A10) 
as as as 
whose roots can be obtained by the familiar formula 
b 1 1 


It is manifest, then, that for each combination of \ and cR, the coefficients ap, a1, @2, and as can be computed from 
Eq. (34), Eq. (A6) can be set up, its roots ;, b2, and b; found, the coefficients N(b,) and gn(b,) evaluated and the 


acceleration response F(s) obtained from Eq. (A4) asa function of s. 


APPENDIX (B).—EFFECcTS OF /7; # 0 AND APPARENT Mass DUE TO WING TwISsT 


In the problem presented here, if the apparent mass due to the wing twist @ and the part of quasi-steady and 
wake forces due to 75 ¥ 0 are considered, if the wing center of twist is located at one-quarter chord length ahead of 
the leading edge of the model wing (i.e., Js) = 1.50 and ls = 2.0), and if the relation Eq. (1) between @(s) and 
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z(s) still holds then the basic equation of motion Eq. (19) becomes as follows: 
MU? (?) 2"(s) = — (*) [W(s) — ¥(0)] pl?Sm — o) + + 
2 
R 


Following the same procedure of solution and using the following operation and initial conditions, in addition to 
those in the text: 


| do — 2"(s) + — + 


= = p%s"(p) — ps"(+0) — | (Be) 
6(+0) = 6’(+0) = 0 
we obtain the response function for this case expressed as (if a5 is not equal to zero) 
F(s = es N(bz) e”” 2S N(6s) els N (6s) 
as qu(be) qu(bs) ga(bs) ds) 
e-0-18 ai 
qe( — 0.130) 
where 4y, be, bs, by, and 0; are five roots of the following fifth degree equation: 
in which 
a; = 0.75cR 
a, = 1.7591cR 
a3 = 1 + + .2446cR . (B5) 
az = 1.3455 + 1.3820 + 0.3149cR 
a, = 0.57525 + 0.0546 + 0.01365cR 
dy = 0.0273 
In Eq. (B3), V(p) is same as that given by Eq. (A2), 
; a4 a3 2 
5 5 
= (p — bi) (b — be) (b — bs) (pb — (p — by (pb + 0.130) (p + 1) (B6) 
Qn(P) = — bn) (B7) 


A sample computation has been carried out for this case for the combination of cR = 75 and \ = 35.60. The 
results of F(s) values as compared with those from the analysis in the text are as follows: 


for s = 0, 5, 10, 20, 40, 60 | 
Appendix (B)F(s) = 0.0000, 0.4899, 0.5552, 0.5677, 0.4861, 0.4012 
Text F(s) 0.0000, 0.5127, 0.5608, 0.5633, 0.4799, 0. 30624 


It is seen that the results based upon these two analyses are so close that, for this particular combination of \ 
and ¢R, it is perfectly justified to use the simplifying assumptions of neglecting the apparent mass due to the 
wing twist and the quasi-steady and wake forces due to /;, ¥ 0. 

Since this combination of \ and cR represents the average case within the ranges of these two parameters investi- 
gated, it seems reasonable to predict that for a trend study these simplifying assumptions are justified for all con- 
ventional aircraft. 


APPENDIX (C).—ESTIMATE OF WING TORSIONAL it is seen that @ may be more generally thought of as the 
FLEXIBILITY mean change of the effective angle of attack due to 

wing elasticity. Then the parameter R may mean to 

In the present analysis, the torsional flexibility of an jnclude the effects on angle of attack of both the tor- 


elastic wing is assumed to be representable by a single sional and the bending flexibility of the actual wing. 
parameter R as defined by Eq. (1). In what follows (Concluded on page 860) 
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Some Stall and Surge Phenomena in Axial- 
Flow Compressors 


MERLE C. HUPPERT* AND WILLIAM A. BENSER?# 
Lewis Flight Propulsion Laboratory, NACA 


SUMMARY 


Observations of rotating stall have shown that a wide variety 
of stall patterns is possible. 

Hot-wire anemometer data on a multistage compressor have 
shown a progressive-type stall at low speeds. The amplitude 
of the flow fluctuations increases in magnitude through the first 
few stages and then diminishes rapidly to a small value in the 
latter stages. A stage-stacking analysis has shown that rotating 
stall will exist over a large portion of the compressor map at low 
speeds but will be instigated almost simultaneously with com- 
pressor surge at high speeds. 

Blades failures attributable to resonant vibrations excited by 
rotating. stall have been experienced in single- and multistage 
compressors. 

In the stage-stacking analysis no deterioration of stage per- 
formance due tounsteady flow resulting from stall of adjacent stages 
was considered. In general, the pressure drop at the stall point 
is believed to be much larger than indicated by an analytical for- 
mulation of compressor performance. Compressor surge is at- 
tributed to a limit cycle operation about the compressor stall 
point, and, as indicated in a few compressor tests and in jet- 
engine tests, a small compressor discharge receiver volume may 
result simply in stall of the compressor without the cyclic char- 
acteristics of compressor surge. In this event, engine operation 
will be limited because of the large drop in performance which 


accompanies compressor stall. 


INTRODUCTION 


i hes GOALS OF INCREASED FUEL ECONOMY and in- 
creased altitude operational ceiling for aircraft jet 
engines have led to engines having compressor pressure 
ratios of the order of 7 to 1 and higher. At the same 
time, minimum engine weight and frontal area are 
essential for operational range of the aircraft. Be- 
cause of its advantages of high flow capacity, high 
efficiency, and simplicity of staging, the axial-flow- 
type compressor is particularly adapted to the high- 
performance jet engine. The stall or surge limit of the 
high-pressure-ratio axial-flow compressor, however, 
presents a major problem with regard to engine accel- 
eration and off-design performance of the engine. This 
limit, which represents the minimum weight flow obtain- 
able consistent with useful operation at any compressor 
rotational speed, is alternately defined as stall limit or 
surge limit. In this paper the stall limit of the multi- 
stage compressor will be referred to as compressor stall 
as differentiated from individual stage stall. 


Presented at the Flight Propulsion Session, Twenty-First 
Annual Meeting, IAS, New York, January 26-29, 1953. 

* Aeronautical Research Scientist. 

+ Chief, Compressor Research Branch. 


The limitations resulting from compressor stall can 
be illustrated by superimposing an engine equilibrium 
operating line on a conventional compressor map (Fig. 
1). In this figure compressor pressure ratio is plotted 
against equivalent weight flow for constant values of 
equivalent or aerodynamic speed. The equilibrium 
operating line shown passes through the region of high 
compressor efficiency at design speed and represents 
steady-state operating conditions for the compressor 
in a jet engine with a fixed exhaust-nozzle area. The 
compressor stall limit is also shown. 

The various off-design performance problems _re- 
sulting from compressor stall can now be discussed. 
The first problem concerns that of accelerating the en- 
gine from idle to full speed. When the fuel flow is in- 
creased from the value required for idling, the turbine- 
inlet temperature increases and results in a transient 
decrease in weight flow. The compressor operating 
point thus moves closer to the compressor stall limit. 
The distance between the equilibrium operating line 
and the compressor stall limit determines the per- 
missible acceleration and is defined as the acceleration 
margin. Large values of this margin permit large in- 
creases in fuel flow and, therefore, rapid engine ac- 
celeration. On the other hand, if the acceleration margin 
is small, large increases in fuel flow result in compressor 
stall or surge, and the engine acceleration is diminished ; 
in fact, in some cases the compressor performance may 
deteriorate so much that the engine will decelerate. 
Stall-limit lines for high-pressure-ratio axial-flow com- 
pressors frequently have kinks that intersect the 
equilibrium operating line and, thus, give zero or nega- 
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Fic. 1. Overall performance of a typical axial-flow compressor. 
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tive acceleration margins at intermediate speeds. 
Twin-spool engines or auxiliary means of acceleration, 
such as adjustable exhaust nozzles, compressor bleed, 
or adjustable compressor-inlet guide vanes that in- 
crease engine weight and complexity, are then re- 
quired in order to obtain satisfactory acceleration. 
The acceleration margin can be increased by overde- 
signing the compressor and turbire, but this technique 
is uneconomical with respect to weight and design- 
speed fuel consumption. The frequently encountered 
adverse intermediate-speed stall characteristics of high- 
pressure-ratio multistage axial-flow compressors thus 
impose a serious acceleration limit for the high- 
performance, high-pressure-ratio jet engine. 

The second off-design problem associated with com- 
pressor stall concerns the development of high engine 
thrust with engines having adjustable jet nozzles or 
other thrust-augmentation devices. Closing the jet 
nozzle and increasing the fuel flow again imposes an 
increased throttling on the compressor, and the oper- 
ating point moves toward the stall limit. The thrust 
increase obtainable is thus determined by the margin 
between the compressor stall-limit line and the cruise 
or open jet nozzle condition. Inasmuch as the peak 
efficiency range for the compressor is usually rela- 
tively close to the stall limit, it is desirable to match 
the engine components so that the compressor operates 
at its peak efficiency for a cruise or open jet nozzle 
condition. The thrust increase by jet nozzle adjust- 
ment may, thus, be seriously limited by the compressor 
stall limit. A similar problem exists for high-thrust, 
high-altitude flight, for which the compressor-inlet 
temperature may be lower than standard sea-level 
temperature. For an engine operating at a constant 
mechanical speed, the aerodynamic or equivalent speed 
is then higher than the design value; and the com- 
pressor stall limit may also be encountered for operation 
of the engine at high thrust, greater than design equiva- 
lent speed. Thus, the compressor stall limit is im- 
portant over most of the engine speed range from idle to 
maximum aerodynamic speed at altitude. 

Because of the importance of the compressor stall 
limit, the source of this limit is a pertinent matter. It 
can readily be shown that stall of the multistage com- 
pressor results from stall of one or more of the individual 
stages; and the stages that instigate overall compressor 
stall will vary with the operating condition, as can be 
seen from the following considerations. At peak effi- 
ciency at design speed, each of the stages of the com- 
pressor will operate at or in the region of its peak per- 
formance. As the pressure ratio is increased at design 
speed, the operating condition of the front stages 
changes only slightly; whereas, the rear stages ap- 
proach their stall points very rapidly because of the 
progressive increase above the design value of density 
of the air passing through the compressor. The small 
change in inlet-stage operating conditions is indicated 
by the nearly vertical characteristic of pressure ratio 
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against weight flow as shown on the compressor map. 
Therefore, at speeds near the design value, stall of the 
multistage compressor is associated with deterioration 
of performance of the latter stages. 

Next the problem of low-speed compressor stall js 
considered. Because of the small stage pressure and 
density ratios associated with low compressor speeds, 
the axial velocity increases excessively through the com- 
pressor, and the rear stages operate at or near their 
choke or turbining condition. The flow restriction of 
the rear stages causes the front stages to operate in a 
high angle of attack or stalled condition. Thus, low- 
speed stall of the multistage compressor is associated 
with deterioration of the performance of the inlet or 
front stages. This is borne out by the fact that the 
flow at 50 per cent of design speed is on the order of 
30 per cent of design speed flow instead of the 50 per 
cent that would be required for design angles of attack 
for the inlet stages. At intermediate speeds, stall will 
be approached more or less simultaneously in several 
stages. It should be noted that increasing the overall 
pressure ratio of the compressor increases the’ devia- 
tions from the design vector diagrams for a given off- 
design speed or weight flow and renders the off-design 
performance and overall compressor stall problems 
more acute. 

Inasmuch as multistage compressor stall is a result 
of individual stage stall, a study of compressor stall 
requires an understanding of individual stage charac- 
teristics and the stage operating conditions for off- 
design operation. For this paper the discussion of this 
phenomenon has been divided into two phases. The 
first is a study of stall characteristics of multistage com- 
pressors and the effects of the individual stage charac- 


teristics on multistage compressor stall. Inasmuch as 


the experimental technique of studying multistage com- 
pressor stall may be hampered by surge and stage 
interactions, it is difficult to evaluate stage character- 
istics in the proximity of overall compressor stall. 
Therefore, much of the investigation reported herein 
has been done by analytical techniques. In addition, 
an analysis of multistage compressor surge charac- 
teristics is considered, and the effects of compressor 
stall and external system characteristics on the compres- 
sor surge phenomenon are discussed. 


STALL OF SINGLE-STAGE AXIAL-FLOW COMPRESSORS 


Single-stage stall characteristics must be evaluated 
in order to understand multistage compressor sta!l char- 
acteristics. Performance of a single-stage axial-flow 
compressor is generally measured with steady-state 
instrumentation that has sufficient damping to attenu- 
ate fluctuations with frequencies over 5 or 10 cycles 
persec. As the flow is decreased beyond the stall point, 
the flow in all blade passages consequently appears to 
be uniform with time. With some stages, regions of 
flow rate exist where surge or low-frequency pulsation 
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is encountered, but steady-state conditions are again 
indicated if further reductions in flow rate are made. 
Recent studies with high-frequency-response instru- 
mentation, such as hot-wire anemometers and pressure 
transducers, have shown stage stall to be a nonsteady 
type of phenomenon in which some sectors of the annu- 
lus are operating at an extremely low mass-flow rate 
as compared with the remainder of the annulus; 
these low-flow sectors rotate about the compressor 
axis. A typical stall pattern is illustrated in Fig. 2 
for a single-stage compressor consisting of guide vanes, 
rotor, and stator. In this example, four stall or low- 
flow zones at the blade tips are depicted by the hazy 
areas. ‘These low-flow zones extend completely through 
the stage and rotate in the direction of the rotor but at 
a lower speed. The mass-flow rate measured on a fixed 
hot-wire anemometer at the discharge of the stage is 
also presented on this figure. This instrument was 
so oriented as to read principally the component of 
flow in the design direction. The alternate high- and 
low-flow zones passing the fixed instrument produce the 
wave form indicated. Similar patterns are observed 
some distance ahead of the guide vanes or behind the 
stators. Temperature and pressure variations across 
the stall zone were also observed but were not evaluated. 

At some point far upstream of the stall zone there 
exists a uniform flow. The stream lines immediately up- 
stream of the blades must, therefore, diverge to main- 
tain continuity around the low-flow or stall zone. This 
spillage of flow to either side of the stall zones, which 
tends to increase the angle of attack on one side of the 
stall and decrease it on the other, results in a tangential 
propagation of the stall zone. The propagation rate 
of the stall region relative to the rotor is generally op- 
posite in direction but lower in magnitude than the 
rotor-blade speed. On the absolute basis, therefore, 
the stall rotates in the direction of the rotor but at a 
lower speed. It appears that in principle the rotating 
or propagating stall may be instigated by a stator as 
well as a rotor. 

To date, no definite relation between the number of 
stall regions in the stall pattern and compressor geom- 
etry has been determined. Two fairly distinct cate- 
gories of stall patterns, however, have been observed: 
one, a progressive stall as shown in Fig. 2, originates at 
the blade ends and consists of multiple stall zones that 
extend over only a portion of the blade span, and the 
other, a root-to-tip stall, generally consists of only one 
zone that extends over the entire blade span. 

Progressive Stage Stall.—The progressive-type stall 
has been observed on several low hub-tip ratio stages 
of the basic solid-body rotation type. For this type 
stage operating with uniform inlet flow, the rotor-tip 
angle of attack varies much more rapidly with flow 
rate than does the root angle of attack, and the rotor- 
blade tip naturally stalls first, while the remainder of 
the blade is inherently unstalled. The circumferential 
and radial extent of the low-flow zones at the blade tip 
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increases with further decreases in flow rate, while the 
root-section operating condition remains relatively un- 
changed. Other types of stage diagram, such as a free- 
vortex flow, might have a progressive-type stall origi- 
nating at the root section of the blade. In general, 
however, the jet-engine requirements of high flow ca- 
pacity and high stage-pressure ratio lead to inlet-stage 
designs in which the rate of change of loading with flow 
at the blade tip is much greater than that at the blade 
root. 

A typical performance curve for a stage that is 
critical at the blade tip is shown in Fig. 3 as a plot of 
pressure coefficient Y against flow coefficient ¢g. As 
stall first occurred (point A), three stall zones were 
found. On the left side of the figure, the oscillograms 
of mass-flow rate obtained by hot-wire anemometers 
for the tip, mean radius, and root section are shown. 
The curves of flow rate against time show the greatest 
flow-rate variations at the tip, a smaller variation at 
the mean radius section, and only minor flow variations 
at the root. As the flow is further reduced from the stall 
point, the area of each stall zone increases. Reductions 
in flow to point B (Fig. 3) resulted in a sudden transition 
from three to four stall zones, and still further reduc- 
tions of flow resulted in a transition to five stall zones 
at point C. Thus, the changes in area of individual stall 
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Fic. 3. Performance of low hub-tip ratio stage with progressive 
stall. 
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Fic. 4. Performance of high hub-tip ratio stage with root-to-tip 
stall. 


zones and changes in number of stall zones result in a 
relatively uniform rate of change of stalled flow area 
with flow rate. Progressive-type stall patterns having 
from three to eight stall zones with rotational speeds of 
40 to 85 per cent of rotor speed have been observed. 
For each stall pattern the rotational speed of the stall 
zones varied directly with rotor speed. It should be 
noted that progressive-type rotating stall results in a 
gradual drop in pressure coefficient as the severity of 
stall increases. 

Root-to-Tip Stage Stall.—When the stall angle of at- 
tack is approached nearly simultaneously over the en- 
tire span of the blade, as in high hub-tip ratio stages, 
the resultant flow-rate variations are relatively uniform 
in intensity from root to tip. Observations of this 
type stall, defined as a root-to-tip stall, show only one 
stall zone. A typical curve of performance in terms 
of pressure coefficient against flow coefficient is shown 
in Fig. 4, together with oscillograms of mass-flow rate 
against time for the tip and root sections. There is 
little radial variation of stall pattern, and from the 
pressure-coefficient curve it is seen that the instigation 
of stall is accompanied by a sharp drop in pressure ratio 
as shown at point A, as compared with the gradual de- 
crease in pressure ratio for the progressive-type stall. 
The drop in pressure rise is on the order of one-third 
of the stage-pressure rise prior to stall. Further reduc- 
tions in flow result in an increase in the circumferential 
extent of the stall zone and a gradual decrease in pres- 
sure coefficient. The stall zone may cover one-third 
or more of the annulus (in circumferential extent) and 
has an absolute rotational speed of 30 to 40 per cent 
of the rotor speed. A root-to-tip stall zone having 
been established, the flow may be considered to be 
divided in two parts, a low-flow or stalled part and a 
high-flow or unstalled part. Changes in compressor 
weight flow may then result in changes in area or weight- 
flow rate of the two zones. In order to unstall this 
type stage, the stalled zone must be completely elimi- 
nated. Because of the freedom of change of area or of 
flow rate in the stalled zones, these zones will persist 
when the compressor weight flow is increased above the 
value for instigation of stall. Unstalling of the stage 
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shown on Fig. 4 is obtained at the point indicated as 
B, and a hysteresis effect results. Thus, the stage 
performance curve has discontinuities between the 
branches. 


The value of the pressure coefficient on the stalled 
branch of the curve in Fig. 4 was determined by con- 
ventional steady-state instruments and represents some 
average value of the pressure coefficient for the stalled 
and unstalled portions of the annulus. Because of the 
response characteristics of the instrumentation used, 
the pressure coefficient in the stall region does not nec- 
essarily represent a true average of the stage perform- 
ance. The existence of this double-branch type of 
stage characteristic is frequently unsuspected or ig- 
nored in the presentation of single-stage performance. 
Discontinuties and double-valued pressure coefficients, 
however, must be considered in computations of multi- 
stage compressor performance from single-stage data 
and are extremely important with regard to studies of 
compressor stall and surge, inasmuch as a double-valued 
single-stage characteristic must effect a double-valued 
multistage characteristic. 


General Discussion of Rotating Stall.—In general, 
low hub-tip ratio stages or stages that approach critical 
conditions much more rapidly at one end of the blade 
span than at the other exhibit progressive-type stall. 
High hub-tip ratio or exit stages in which angle of 
attack variations along the blade span are small exhibit 
root-to-tip stall. Available data on _ intermediate 
hub-tip ratio or typical middle stages of a multistage 
compressor exhibit first a progressive stall that, with 
further reductions in flow, changes to a root-to-tip 
stall with the associated sharp drop in pressure and 
the hysteresis characteristic. An important difference 
between progressive and root-to-tip rotating stall is 
the gradual decrease in pressure coefficient associated 
with progressive stall as compared with the discon- 
tinuity of pressure coefficient associated with root-to- 
tip stall. Of secondary importance is the fact that the 
frequencies of flow variation are generally much higher 
for progressive than for root-to-tip stall. 


The rotational velocities of the stall zones are ap- 
proximately proportional to the compressor speed, and 
in the useful speed range the multiple-zone progressive 
stall may result in flow-rate variations with frequencies 
relative to the rotor or stator blading which are of the 
same order as the resonant bending frequencies of 
typical inlet-stage blading. Inasmuch as frequent 
blade failures have resulted from brief periods of oper- 
ation in stall at intermediate compressor speeds, indi- 
cations are that rotating stall is a possible source of 
resonant excitation that may lead to blade failure. 


It should be noted that the rotating-stall pattern is 
not always as sharply defined as indicated in the stages 
discussed, and at least one stage has been tested in which 
no rotating stall existed. This particular stage has an 
extremely high camber at the blade tip. Thus, although 
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STALL AND SURGE 
rotating stall appears to be prevalent, it is not neces- 


sarily the only type of stall which can exist. 


STALL OF MULTISTAGE COMPRESSORS 


The second phase of this study of axial-flow- 
compressor stall is consideration of the effects of the 
single-stage stall characteristics on the multistage stall 
problem. As was discussed earlier in this paper, for 
operation of high-pressure-ratio multistage compressors 
at 50 per cent of design speed, flow restrictions in the 
latter stages limit the flow and result in stalling of the 
front stages. In spite of front-stage stalling, however, 
the compressors do have an appreciable flow range of 
usable operation at low speeds. It is evident, there- 
fore, that stage stall is a necessary but not a sufficient 
condition for multistage compressor stall. 

Experimental Multistage Study.._Some understanding 
of the flow processes prevalent at low and intermediate 
speeds is given by some limited hot-wire anemometer 
data recently taken on a multistage research com- 
pressor. These data indicate that a rotating stall of 
the progressive type existed during the surge-free and 
useful range of compressor operation at low speeds. 
A summary of the radial and axial extent of these ro- 
tating-stall patterns is given in Fig. 5. The ordinate 
of the curves is the ratio of the amplitude of the flow 
variations Ap!’ of rotating stall to the average pV for 
the tip, mean, and root sections. The abscissa is the 
axial distance through the compressor. A _ rotating 
stall of the progressive type existed in the earlier stages. 
The amplitude of the flow variations in the stall zones 
increased through the first few stages and then de- 
creased rapidly to relatively small values in the rear 
stages. The amplitude of the flow variations de- 
creases rapidly with radius and is relatively small at the 
root section even in the inlet stages. Stall patterns 
having three to seven stalled zones were found with 
the number of stalls increasing with a decrease in flow 
and speed. A marked similarity exists between the 
stall patterns measured in the front stages of this com- 
pressor at low speed and those obtained with the inlet- 
type single-stage compressor discussed previously. 
The stall zones extended nearly axially through the 
compressor. and the harmonics of the flow variations 
were damped in the latter stages, leaving a sinusoidal 
pattern at the compressor discharge. As was men- 
tioned: in the: single-stage discussion, the stall frequen- 
cies measured for this progressive-type stall were of 
the same order of magnitude as the natural bending 
frequency of the blading of the first few stages. There- 
fore, resonant excitation of rotor and stator blades in 
the front stages as a result of rotating stall is a potential 
source of blade failure. Although these studies were of 
a preliminary nature and the results are only qualita- 
tive, they clearly indicate that, even though rotating 
Stalls instigated by the inlet stages persist throughout 
the compressor, the pressure-ratio producing capacity 
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of the compressor is not excessively inhibited so long 
as the flow coefficients of middle and latter stages are 
sufficiently above their stall values. 

Computed Multistage Performance.—The results of 
these preliminary experimental studies suggest an 
analytical study of overall multistage compressor stall 
and off-design operation based on a knowledge of indi- 
vidual stage characteristics. An approximate knowl- 
edge of variations in operating conditions of each 
stage can be obtained over the complete range of multi- 
stage operation by a stage-stacking technique. A 
representative analysis was made, therefore, for a 
hypothetical 12-stage compressor having an overall 
pressure ratio: of 7.75 at the reference or match point. 
For this analysis, stage performance was generalized 
in the form of pressure coefficient and adiabatic effi- 
ciency against a flow coefficient based on mean radius 
velocities, and the overall compressor performance was 
computed with a stage-by-stage technique. Single- 
stage data suggest that high-flow, high-pressure-ratio 
inlet stages will probably have progressive stall; whereas 
middle stages will have progressive stall followed by 
root-to-tip stall, and exit stages will have only root- 
to-tip stall. The first four stages, therefore, were 
assumed to have a progressive stall, and the stage 
characteristics shown in Fig. 6 were chosen. For 
stages 5 through 8, an initial progressive stall followed 
by a root-to-tip stall at a lower flow coefficient was 
assumed; this assumption gave the stage character- 
istics indicated. The last four stages were assumed to 
have only a root-to-tip stall, and thus the selection of 
the rear-stage-performance characteristics was as 
shown. 

The curves for all stages were identical except for the 
stall characteristics. A reference speed was chosen and 
the area ratio through the compressor determined so 
that all stages operated at a flow coefficient of 0.69 at the 
reference condition. Overall compressor performance 
was then computed for speeds from 50 to 100 per cent 


1 as 
the 
lled 
‘on- 
lled 
the | 
ed, 
‘m- 
of 
ig- 
ce, 
Its, 
Iti- 
ita 
of 
ed 
ed 
al, 
al 
de 
Il. 
of 
vit 
te 
xe 
th 
p 
d 
iS 
d 
e 
| 


> 
1 


w 


JOURNAL OF THE AERONAUTICAL 

staces 1-4 [ STAGES 5-8 [ STAGES 9-12 
ae a 
2 
a 
= & 


PRESSURE COEFFICIENT, 


FLOW COEFFICIENT, NACA 


Stage performance parameters for stacking analysis. 


PRESSURE RATIO, Pi3/P| 
D 
T i T 


T 


a 
80 “PERCENT OF 
/ REFERENCE 
/ 80 SPEED 


a 8 12 16 20 24 28 32 36 


SPECIFIC WEIGHT FLOW, W/8/8A 


Fic. 7. Computed multistage compressor performance. 
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of the reference speed in 10 per cent intervals. It 
should be noted that the stacking technique based on 
generalized stage curves does not consider decreases in 
either the stage range or efficiency at extremely high 
Mach Numbers, nor does it account for changes in 
performance of a given stage as a result of nonsteady or 
nonuniform flow imposed by adjacent stages. This 
technique is used primarily to determine the relative 
matching of the stages for off-design performance points 
for the muitistage compressor and, thus, to determine 
how peculiarities in stage characteristics resulting from 
instigation of rotating stall may affect overall perform- 
ance. 

The computed compressor map in terms of pressure 
ratio against specific weight flow (flow rate per unit 
of annulus area) is presented in Fig. 7. An appre- 
ciable difference in peak efficiency at high and low 
speeds is indicated by the contours of constant effi- 
ciency. Efficiencies of slightly over 0.70 were com- 
puted for 50 per cent speed, whereas a peak efficiency 
of over 0.85 was obtained for the reference speed. The 
islands of high efficiency at the reference speed are 
broader than those generally measured experimentally; 
this difference can be attributed to the fact that neither 
the range of high individual stage efficiency nor thie 
maximum of stage efficiency was reduced as Mach Num- 
ber was increased. 

The lowest flow rate considered at each speed was 
the one at which some stage had a discontinuity in 
pressure coefficient resulting from the instigation of 
root-to-tip stall. For these minimum flow points, 
the overall pressure ratio of the multistage compressor 
was computed for both the upper and lower branches 
of the individual-stage curve. For each speed, there- 
fore, a discontinuity was computed for the overall 
multistage pressure-ratio curve, as shown in Fig. 8. 

Stagewise Variation of Flow Coeffictent.—One of the 
principal results of the stage-stacking analysis is the 
determination of the particular stage responsible for 
the discontinuity in overall pressure ratio at each speed. 
The study of the sources of these discontinuities is pro- 
vided by plotting the stagewise variation of flow co- 
efficient for the points where the discontinuities occur 
for three representative speeds, 50, SO, and 100 per cent 
of the reference value (Fig. 8). The originally as- 
sumed pressure-coefficient curves are superimposed at 
the top of the figure, and the reference flow coefficient 
of 0.69 for all stages is indicated. The progressive 
stall limit for stages 1 through 8, at a flow coefficient of 
0.565, is indicated by the dashed line on Fig. 8. The 
root-to-tip stall for stages 5 through 8, at a value of 
flow coefficient of 0.470, and the associated discon- 
tinuity in pressure coefficient are indicated in the curve 
at the top of the figure. Similarly, the discontinuity 
of pressure coefficient resulting from root-to-tip stall 
for stages 9 through 12 occurs at a value of flow co- 
efficient of 0.565. The flow coefficient at the inlet 
stage for the reference-speed discontinuity point is only 
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STALL AND SURGE 
slightly lower than the value used in the reference point 
computation. Succeeding stages, however, have a pro- 
gressive decrease in flow coefficient from the value for 
the inlet stage, and the stall of the twelfth stage is re- 
sponsible for the discontinuity in the overall pressure- 
ratio curve that was shown on the previous figure. 
Effects of compressibility and radial variations on stage 
range and efficiency, which were neglected in this anal- 
ysis, may reduce the maximum pressure ratio from 
the value computed here, but the general shape of the 
overall pressure-ratio curve and the determination of 
the instigating stage are valid. 

As compressor speed is reduced, there is a general 
decrease in flow coefficient for the inlet stages and an 
increase for the rear stages. This is borne out by the 
stagewise variations of flow coefficient for the discon- 
tinuity point at SO per cent of reference speed. At 
this speed all stages up to stage 9 are operating in the 
range where progressive rotating stall without discon- 
tinuities exists, and the discontinuity in the overall 
performance curve results from instigation of root-to- 
tip stall in the ninth stage. Again, consideration of 
effects of root-to-tip stall on the performance of other 
stages would result in a large drop in the attainable 
overall pressure ratio at SO per cent of reference speed. 

For the 50 per cent speed point, the trend is similar 
to that for SO per cent speed, except that the discon- 
tinuity of pressure rise results from the instigation of 
root-to-tip stall in the fifth stage. As a result of the 
progressive stagewise variations of flow coefficient and 
the grouping of stage types, root-to-tip stall was en- 
countered only in stages 5, 9, and 12. With other 
designs, this situation will, of course, be altered. In 
general, however, early stages will instigate stall at 
low speeds, and latter stages at high speeds. 

In order to study the speed ranges over which various 
stages are stalled, the overall performance map has been 
replotted (Fig. 9), and lines of first-, fifth-, ninth-, and 
twelfth-stage stall have been added. The discon- 
tinuities in overall pressure ratio for speeds of 50, 60, 
and 70 per cent of the reference value are the result of 
instigation of root-to-tip stall in the fifth stage. For 
speeds of SO and 90 per cent, the discontinuity is caused 
by the ninth stage. At higher speeds, root-to-tip stall 
of the twelfth stage causes the discontinuity in the 
overall pressure-ratio curve. The magnitude of pres- 
sure drop shown is the result of stall of only one stage. 
As a result of flow fluctuations resulting from root-to- 
tip stall, however, adjacent stages that are close to 
stall would probably be adversely affected, and the 
magnitude of pressure discontinuity would be much 
greater than indicated. 

All points to the left of the inlet-stage stall-limit line 
indicate the existence of a progressive stall in one or 
more of the inlet stages. All the flow range for 50, 60, 
and 70 per cent speeds and most of the flow range for 
80 per cent speed are in the range of inlet-stage stall. 
The fifth-stage progressive-stall limit line shows that 
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Fic. 9. Stage stall limits for multistage compressor. 
progressive stall exists at least back to the fifth stage 
over the low-flow range for speed up to 90 per cent of 
reference. As shown on Fig. 9, therefore, progressive 
stall exists in many of the inlet stages at low speeds for 
all flows, and it also exists for low flows at intermediate 
speeds. Discontinuities in overall pressure are in- 
stigated by root-to-tip stall of the middle stages at low 
speed and the latter stages at high speed. 

General Considerations.Effects of progressive-type 
stall on altering the effective performance curves of 
succeeding stages have not been considered in this 
analysis. Such interactions would probably result in 
some deterioration of stage performance but would 
not appreciably alter the general trends of the pressure 
ratio against weight flow at a constant speed. At in- 
termediate speeds, where a large number of stages are 
approaching stall simultaneously as flow is decreased, 
root-to-tip stall of middle stages may occur at a higher 
value of flow coefficient than normally as a result of the 
progressive stall of early stages. This situation would 
cause the discontinuity of overall pressure ratio to oc- 
cur at a higher value of specific weight flow and at a 
lower pressure ratio than indicated by this analysis. 

The relations of the sharp discontinuities of pres- 
sure ratio shown to compressor surge on Fig. 9 will be 
discussed later; but whether or not the compressor 
surges, these points represent a lower limit to the useful 
flow range at any speed. It will be noted that, at speeds 
of 80 and 90 per cent, the discontinuities in pressure 
ratio occur at values of pressure ratio less than the peak. 
On the basis of stability criteria, surge may be encoun- 
tered at approximately the flow for peak pressure ratio; 
therefore, a curve through the peak pressure-ratio 
points at various speeds has been included in Fig. 9. 
The curves of neither overall stall point nor peak pres- 
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sure ratio indicate the kinks in the compressor stall or 
surge line which are prevalent for high-pressure-ratio 
axial-flow compressors. Observation of regions where 
progressive-type rotating stall prevails, however, sug- 
gests that stage interactions may induce root-to-tip 
stall in inlet stages, which might readily cause kinks 
in the compressor stall-limit line. Variations in stage 
matching and stage characteristics may also be re- 
sponsible for the undesirable kinks. As indicated in 
the analysis presented herein, a large group of stages 
in the middle of the compressor are operating near their 
individual stall points at intermediate speeds; and, 
if any one stage stalls, a serious deterioration of many 
stages may be precipitated. Further studies of stages 
interactions are necessary to predict completely com- 
pressor stall characteristics at low and intermediate 
speeds. A stall may cause discontinuities in indi- 
vidual stage characteristics and result in sharp dis- 
continuities in the overall compressor performance 
curve. A definite limit of useful operating range of 
the compressor at any speed, regardless of whether 
surge occurs, is thus directly attributable to the dis- 
continuity in the performance of one or more of the 
individual stages. 


COMPRESSOR SURGE 


Tests of multistage axial-flow compressors are con- 
ventionally made in test rigs having large inlet and dis- 
charge tanks. Violent surge is commonly encountered 
in these rigs as the compressor weight flow is decreased 
to a limiting value. This surge is a severe audible 
pulsation of compressor pressure rise and weight-flow 
rate and normally has a frequency on the order of 0.5 
to 2 cycles per sec., although much higher values have 
been observed on small rigs. In some jet-engine tests 
and in multistage tests with low receiver volumes, stall 
with a severe drop in pressure ratio and efficiency has 
been observed without the audible flow pulsation nor- 
mally defined as surge. Many investigators have 
studied the criteria for surge as the criteria for sta- 
bility, considering the response of the system to small 
variations from an equilibrium operating point. A con- 
tinuous pressure-ratio, weigit-flow characteristic was 
implied in these analyses." ? 

Classical Surge.—In the previously reported investi- 
gations, in which the compressor characteristic of pres- 
sure rise was assumed to be a continuous function of 
weight flow for a constant compressor speed, lineariza- 
tion of the compressor and discharge throttle character- 
istics and simplification of the external system result in a 
solution of the differential equation of flow of the fol- 
lowing form: 


AW = sin (wt + (1) 


where 


AW = the variation of weight flow from the equilib- 
rium value 
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Ki = aconustant 

a = damping coefficient 

t = time 

= a function of the frequency of oscillation 
AK, = aconstant expressing the phase relation 


From this equation it can readily be seen that the 
damping coefficient @ defines the stability of the sys- 
tem— that is, its response to a small deviation from the 
equilibrium operating condition. If the damping 
coefficient has a positive value, the system will be 
stable. If the damping coefficient has a negative value, 
however, the system will be unstable and the magnitude 
of the resulting oscillations will continue to increase 
until the original assumptions of linearity are no longer 
valid. The effects of the neglected nonlinearity, 
therefore, will limit the amplitude of oscillation to 
some finite value, and the cyclic pulsation of surge re- 
sults. For convenience, this type surge will be referred 
to as classical in this paper. It should be noted that 
the analyses of references 1 and 2 merely define the 
conditions of stability. Determination of amplitudes 
of flow pulsation are dependent on nonlinearities of the 
compressor characteristics in the stall region, and the 
procurement of accurate data of this type is difficult 
because of the existence of compressor surge and ro- 
tating stall. When compressor characteristics are dis- 
continuous, however, a different approach to the anal- 
ysis of surge may be utilized. 

Limit-Cycle Surge.—As discussed earlier this 
paper, single-stage tests have shown pressure-ratio, 
weight-flow characteristics that are discontinuous at 
stall (Fig. 4). This characteristic indicates that, as 
weight flow is even slightly decreased below the value 
for stall, the pressure ratio drops discontinuously to the 
lower branch of the curve. Further decreases in 
weight flow, however, cause a continuous change in 
pressure ratio. Once stall has been encountered, the 
flow must be increased beyond the initial value for 
stall in order to relieve the stalled condition. This re- 
sults in the hysteresis effect shown. Discontinuities 
in stage curves cause discontinuities in overall pressure- 
ratio curves for multistage compressors, as indicated 
in the calculated performance map presented. Double- 
branch performance curves have been obtained in one 
multistage compressor test in which a small discharge 
receiver volume was used and also in a research jet 
engine.’ Operation on the stalled branch of the com- 
pressor characteristic apparently did not result in the 
classical instabilities of surge. These compressors, 
however, surged violently when operated in standard 
component test rigs that had large receiver volumes. 
Inasmuch as the compressor characteristics on the 
stalled branch of the curve apparently satisfied the sta- 
bility requirements of Eq. (1), the question arises as to 
the reason for violent surge when large receiver volumes 
were utilized. The discontinuity in pressure associated 
with stall will undoubtedly result in a damped, transient 
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STALL AND SURGE 


flow-rate fluctuation if normal stability requirements 
are satisfied. If, however, the compressor flow rate 
reaches or exceeds the value for stall recovery during 
the transient, the stall condition will be relieved, and 
the compressor will again seek equilibrium at the orig- 
inal point of stall on the upper or unstalled branch 
of its characteristic. Cyclic stalling and unstalling 
of the compressor, resulting from transient flow vari- 
ations following an initial stall, thus provide a mecha- 
nism for surge. The discontinuous-type compressor 
characteristic, therefore, suggests solving the differential 
equation of motion for the compressor and external 
system for the transient flow-rate response to a step 
change in compressor discharge pressure at stall. 

For the analysis of surge in a compressor with a dis- 
continuous pressure-ratio, weight-flow characteristic, 
the simplified system indicated in Fig. 10 is considered 
in which there are a compressor, discharge receiver, and 
discharge throttle. The compressor draws air from 
the atmosphere at station 1 and discharges to the re- 
ceiver, Station 2. The flow is then discharged through 
the outlet throttle to station 3. Similar to the tech- 
nique of references 1 and 2, the compressor character- 
istics and throttle characteristics are linearized as shown 
in Fig. 11. The upper or unstalled branch of the pres- 
sure characteristic has a negative slope and stalls at 
point 0. The pressure drop at stall is represented by 
AP,, and the lower or stalled branch of the character- 
istic curve is assumed to have a zero slope. For this 
analysis, the throttle characteristic passes through the 
compressor stall point indicated by 0 in Fig. 11. The 
increment of weight flow between the stall point and 
the intersection of the throttle line and the stalled 
branch of the compressor curve is defined as AIV,. 
The flow must be increased above the original stall 
value by an increment AlV,, to permit stall recovery. 
On the basis of the surge criteria previously discussed, 
this system is stable on the upper or unstalled branch 
of the compressor characteristic except at the stall 
point 0. The system is also stable on the lower or 
stalled branch of the compressor characteristic. In 
the solution of the system differential equation, the com- 
pressor is assumed to be operating at the stall point, point 
0 on Fig. 11, and a step change in pressure ratio is im- 
posed. This results in a receiver pressure that is greater 
than the compressor is capable of producing in the 
stalled condition, and the solution of the linearized 
flow equations for the resulting flow variation with time 
is of the following form: 


AW = — AIV. [1 + e~™ (Asin wt — cos wf)} (2) 
where 
AlV == the variation of weight flow from the initial 


or stall value 
absolute value of the weight-flow incre- 
ment between equilibrium points on the 


AW, 
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Fic. 11. Linearized compressor and throttle characteristics 
< 
> 
oO 
] 
> 
° 
0 2 4 6 
TIME, Qt 
Fic. 12. Transient flow-rate fluctuation following stall. 
upper and lower branches of the com- 
pressor characteristic 
a = the damping coefficient determined by 
compressor and system characteristics 
t = time 
A = a constant determined by compressor 
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w = a system constant that determines the 
frequency of oscillations 


It may be noted that if the weight flow at the inter- 
section of the throttle curve and the stalled branch of 
the compressor curve is taken as the base for deter- 
mination of All, the above equation may be expressed 
in the same form as Eq. 1. As in references 1 and 2, 
the magnitude and stability of the resulting flow oscil- 
lation are dependent on the magnitude and sign of the 
damping coefficient a. 

Solutions to the above equation have been computed 
for various compressor and system characteristics, 
and the results are presented with dimensionless vari- 
ables in Fig. 12. The ordinate for this plot is AV + 
|AW,|, or the ratio of the weight-flow increment at any 
instant measured from the equilibrium point on the un- 
stalled branch of the compressor curve to the absolute 
value of the difference in weight flow between the 
equilibrium point on the unstalled and the stalled 
branch. The variable of time is made dimensionless 
by multiplying it by a characteristic frequency (or 
dividing it by a characteristic time). The undamped 
natural frequency 2 of the compressor and _ receiver 
acting as a resonator was selected as the characteristic 
frequency for this presentation, and, accordingly, Qt is 
used as the abscissa for Fig. 12. Similarly, the damp- 
ing term a, which has the dimensions of frequency and 
is inversely proportional to the time constant of the re- 
ceiver and throttle, is divided by 2 to provide the di- 
mensionless parameter a/Q. The variation of AIV + 
|AW,| with Qt is presented for values of a/Q of 0.1, 0.3, 
0.5, and 0.7 in Fig. 12. Negative values of the ratio 
AW/AIV, represent decreases in flow from the initial 
or stall point, and positive values represent increases. 
For small damping, a/Q equal to 0.1, the flow decreases, 
then increases to a value greater than the stall value, 
and then decreases again. If the flow increases above 
the value for stall recovery AIV,,/ IAW, I, the com- 
pressor’s ability to develop pressure will increase dis- 
continuously to the upper or unstalled branch of the 
curve. A similar computation for a step increase in 
pressure would indicate a sharp increase in flow rate, 
followed by a decrease to the value for stall, and the 
cycle would repeat. Inasmuch as the main purpose of 
this computation is to study the potential of surge due 
to cyclic stalling and unstalling, only the overshoot on 
the stalled branch of the compressor curve was con- 
sidered, and the transient following the step increase 
in pressure resulting from stall recovery was not com- 
puted. 

This cyclic-type operating condition, caused by a 
step change in pressure ratio due to complete com- 
pressor stall, has been defined as limit-cycle surge. 
This type surge will result in violent oscillations of flow 
and pressure ratio. Because of inertia and capacitance 
effects of the system, the receiver pressure change will 
lag the flow change during the transient. By again re- 


ferring to Fig. 12, it is seen that, for larger values of 
a/Q such as 0.3 to 0.7 as shown on this chart, the 
transient flow fluctuations damp out rapidly and the 
compressor operating point converges on the intersection 
of the throttle line and the stalled branch of the compres- 
sor characteristic (AW/'AW,| = — 1), and limit-cycle surge 
will not result. From a consideration of discontinuous 
pressure characteristics, therefore, limit-cycle surge 
due to complete compressor stall may or may not occur, 
depending on the magnitude of the system damping. 

It can be seen that two types of surge are possible: 
classical surge resulting from system instabilities for 
a compressor with a continuous characteristic curve, 
and limit-cycle surge resulting from a discontinuity 
of the compressor characteristic. If the system damp- 
ing is sufficiently high at the equilibrium point on the 
stalled branch, compressors with discontinuous char- 
acteristics may simply stall without surging. Classical 
surge may, however, exist as a result of negative damp- 
ing on either the unstalled or stalled branch of a double- 
valued compressor characteristic, and the resulting 
instabilities may incite a limit-cycle-type surge. 

Effects of System Characteristics on Damping Coeffi- 
cient.—From the preceding analysis, the importance 
of system damping on compressor surge or stall char- 
acteristics is evident. Consideration must, there- 
fore, be made of the effects of system characteristics 
on the magnitude of the damping coefficient. If the 
receiver and throttle have a large time constant—that 
is, if the discharge volume is large or the throttle curve 
slope is large, or both, as is the case for most compressor 
component test rigs—the damping coefficient will be 
small, and the compressor will probably surge. It 
should also be noted that limit-cycle surge may be 
encountered more readily at high compressor speeds 
than at low speeds, because of the increase in slope of 
the throttle curve and also because of the very narrow 
range of weight flows for which the two branches over- 
lap. Classical-type surge is also more likely to occur 
at high than at low speeds. 

In jet engines the discharge receiver volume is the 
combustion-chamber volume and is normally small 
when compared to the receivers used in compressor 
test equipment. Moreover, the throttling is done 
by the turbine nozzle, and the throttle curve is essen- 
tially a constant fuel-flow line which has somewhat less 
slope than the throttles used in compressor test facili- 
ties. As a consequence, the compressor in the jet en- 
gine will have a relatively large damping coefficient and 
will tend to stall rather than to surge. 

Experimental Considerations. Further studies of 
multistage compressor characteristics at the stall point 
and in the stall region are vital for obtaining a better 
understanding of surge and stall phenomena. Cau- 
tion, however, must be exercised in obtaining and an- 
alyzing these data. During operation wit overall 
compressor stall, fast-response instruments at either 
the compressor inlet or discharge may indicate large 
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flow variations due to rotating stall, and these fluctu- 
ations may easily be mistaken for compressor surge. 
Phase relations from two fast-response instruments in 
a plane normal to the axis of rotation of the compressor 
will readily identify rotating stall. Furtherinore, fluc- 
tuations in torque required to drive the compressor will 
exist during surge, whereas the torque will be relatively 
constant with stall. 

Inasmuch as the flow fluctuations of surge prevent 
the experimental determination of the slope of the 
compressor characteristic curve at the stall point, a 
discontinuity in the compressor characteristic can be 
found only if surge is prevented— that is, the stall char- 
acteristics of a compressor cannot be properly deter- 
mined if stall results in surge. 


CONCLUSIONS 


From the preceding discussion of the various aspects 
of stall and surge in axial-flow compressors, the following 
results are obtained: Stall of single-stage co:npressors 
frequently results in zones of high and low flow which 
rotate about the compressor axis. These rotating stalls 
may be classed as either progressive stall with multiple- 
stall zones, or root-to-tip stall, which generally con- 
sists of a single-stall zone. The progressive type re- 
sults in a gradual decrease in pressure ratio with flow 
decreases below the stall value and is most likely to 
occur in low hub-tip ratio stages where stalling weight 
flow is significantly greater for one part of the blade 
span than the other. Root-to-tip stall results in sharp 
drops or discontinuities in the pressure-ratio flow char- 
acteristic and generally exists in stages of high hub-tip 
ratio where critical conditions are approached simnul- 
taneously over the entire span of the blade. Rotating 
stall of either type may excite compressor blading to 
resonant vibration. 

Rotating stall has been observed in multistage, as 
well as single-stage, axial-flow compressors. A stacking 
analysis for a hypothetical compressor with front stages 
having a progressive-type stall and with middle and 
latter stages having the pressure-ratio discontinuities 
associated with root-to-tip stall indicates multistage 
compressor performance with discontinuities in the 
pressure-ratio, weight-flow characteristics at all speeds. 
At low speeds these discontinuities result from stall 
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of the earliest stages having discontinuous characteris- 
tics, and at high speeds from stall of the latter stages. 
Progressive-type stall of the inlet stages may exist over 
the entire flow range at low speeds. At intermediate 
speeds all stages are close to their stall points, and 
stage interactions resulting from the flow variations of 
stall of any element may seriously limit the pressure- 
ratio-producing capacity of the compressor. 

In addition to surge that results from classical system 
instabilities, a limit-cycle-type surge may result from 
discontinuities in the pressure-ratio, flow character- 
istic of the compressor. A compressor may stall rather 
than surge, provided the damping of the compressor 
and external system is sufficiently high. As was the 
case for classical surge, the potential of limit-cycle 
surge is greatest at high compressor rotational speeds 
and pressure ratios. The potential of limit-cycle surge 
is also greatest where the slope of the throttle char- 
acteristic is large, the incremental flow increase for 
stall recovery is small, and the discharge receiver vol- 
ume is large, as in compressor component test rigs. 
In jet engines the receiver volume is small, the slope 
of the throttle characteristic is small, and the pressure 
ratio at which stall will normally be encountered is 
low;. consequently, compressors in jet engines are less 
likely to surge than are those operated in test rigs. 

This paper deals primarily with the definition of the 
problem of stall and surge. Much research must still 
be done on such problems as correlations of types of 
stall with stage geometry, evaluations of stage inter- 
actions when rotating stall exists in multistage com- 
pressors, and evaluation of the system characteristics 


governing surge. 
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RIEF REPORTS of investigations in the aeronautical sciences and discussions of papers pub- 
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A Note on Organ-Pipe Resonance in Ducted 
Burners 


J.C. Truman* and N. J. Lipstein 
General Engineering Laboratory, General Electric Company 
Schenectady, N.Y. 


August 21, 1953 


acini OR ORGAN-PIPE VIBRATIONS of the gases in 
ducted burners (such as rocket, ram-jet, and turbojet com- 
bustion systems) have been treated in several recent papers.'~° 
In such systems, the axial velocity of the gas through the burner 
is subsonic but is not necessarily small compared with the ve- 
locity of sound. To predict the resonant frequencies of a given 
system, therefore, it is necessary to investigate the effect of axial 
velocity on the formation of longitudinal standing waves. This 
topic is considered by treating the case of one-dimensional flow in 
a constant diameter duct. 
The following notation will be used: 


e = velocity of sound relative to the medium 
Ci, = constants - 

C3, Cs 
Ei = frequency, cycles per sec. 
fr = resonant frequency for (# = 0), cycles per sec. 
ys = length of duct 
M = Mach Number = 1/c 
n 
p = pressure 

= time 

u = velocity 
x = distance coordinate along duct axis 
p = density 
0 = 1/(1 — M?*) 


Subscript zero indicates steady-state quantities. 

Superscript prime indicates small perturbation quantities 
associated with a plane pressure wave propagated along the duct 
axis. 

First consider the case of zero axial velocity of the gas. The 
instantaneous pressure, density, end velocity can be written in 
terms of steady-state and perturbation components. 


Pp = po + p’ a (1) 
u= 1’ \ 
Applying Eqs. (1) to the equation of continuity 
(Op/dt) + (0/dx) (pu) = 0 (2) 


and Euler’s equation 


* Now Assistant Professor of Aeronautical Engineering, USAF Institute 
of Technology, Wright-Patterson Air Force Base, Ohio. 


Ou ou + Op 
p pu =0 (3) 
ot Ox ot 
and neglecting all terms of second or higher order in the perturba- 
tion quantities results in the wave equation 


(07p'/ot?) — = 0 (4) 


A solution of Eq. (4) may be written 
pb’ = C,] sin — (x — ct) + sin (x + ct) ] + 
2af 
C, | cos — (x — ct) + cos — (x + ct) (5) 


As boundary conditions, let a forcing function of the form 
pb’ = P’ cos 


be maintained at one end (x = 0) of the duct, and let the condi- 
tion (p’ = 0), approximating an open end, exist at (x = L),. 
Rearranging Eq. (5) and applying the assumed boundary condi- 
tions results in 

P* 


sin — (L — x) cos 2aft (6) 


, 


sin (22fL./c) 
This represents a standing wave of amplitude P’/sin (2xfL/c). 
Resonance occurs, therefore, when (22fL./c = nw)—that is, the 

resonant frequency is 
fr = nc/2L (7) 


The above theory is well known, and serves as an approach to 
the case in which the gas flows through the duct with a uniform 
axial velocity uw <c. Eqs. (1) now become 


p Po + 
p= pot p’ (8) 
u = + 


and the corresponding ‘‘wave equation” is 


p’ Op’ 
id + 2% P — (c? — U2) (9) 
of? Oxot Ox? 


A solution of Eq. (9) is 


€ 


C+ 
[x + — ged [x — (¢ + wo)t] + 
cos———_ [xn + (¢c — (10) 
c= the 


Applying the same boundary conditions as before to Eq. (10) 
and introducing M = u/c, we find 
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sin 


cos — x + 2zxft (11) 
é 


where 6 = 1/(1 — V?). 

Comparing Eq. (11) with Eq. (6), it is seen that the condition 
now existing is no longer that of a true standing wave, since a 
traveling component appears in the cosine term. From the 
amplitude term of Eq. (11), resonance is seen to occur when 
(2xfoL/c = nm). The resonant frequency is therefore 


f = nc/2L0 = fr(1 — M?) (12) 


Eq. (12) shows that the effect of an axial velocity through the 
duct is to reduce the resonant organ pipe frequency by a factor 
of (1 — 2). This result was also derived by Trimmer’ using an 
analysis based on the acoustic impedance of pipes, and it was 
stated by Fenn, Forney, and Garmon.? 

The effect of an axial velocity through the duct on the end cor- 
rection associated with an open end pipe has not been investi- 
gated. However, one might expect the usual correction factor 
to be increased as a result of the flow out of the duct. Thus, an 
effectively longer duct would be established, and correspondingly 
a still further decrease in resonant frequency would result. 
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A Note on the Effect of Diffusion Fields on the 
Laminar Boundary Layer 


John W. Smith 
Aerodynamicist, Bell Aircraft Corporation, Buffalo, N.Y. 
August 3, 1953 


ier CHARACTERISTICS OF A BOUNDARY LAYER are strongly 
affected by an admixture of foreign gases, provided their 
physical properties differ markedly from those of the primary 
gas. A theory has been developed! which describes the effect 
of a general diffusion field on the dynamic and thermal char- 
acteristics of a laminar boundary layer on a flat plate in steady 
compressible flow. These characteristics are still determined 
by the well-known boundary-layer equations, 


Ou Ou ou 
pu + pv = (. =) (1) 
op/oy = 0 (2) 
(0/Ox) (pu) + (0/dy) (pv) = 0 (3) 
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where, however, the quantities p (density), u (coefficient of vis- 
cosity), & (thermal conductivity), and c, (specific heat) are now 
functions of the local concentration. For this reason it becomes 
necessary to introduce an additional equation describing the 
diffusion field, before the given system can be integrated. 

The diffusion equation expressing the conservation law for a 
foreign gas can be written in the general form 

on, on, 
v 


= —div J, (5) 
ox oy 


1 


where ”, represents the local particle density of foreign gas and 
J, represents the corresponding diffusion current. In the ab- 


sence of body forces and pressure gradients, one obtains” * the 


following expression for J;: 


— 


J, = —nD,.| grad myo — - r grad 7 (6) 


Here represents the total particle density; Dj», the coeflicient of 
mutual diffusion; m9, the ratio of m to n; and a, the thermal 
diffusion constant. The first term of Eq. (6) represents a diffusion 
current arising from the concentration gradient, while the second 
or thermal diffusion term expresses the tendency of lighter 
molecules to migrate into regions of high temperature. The final 
form of the diffusion equation is obtained by combining Eqs. (5) 
and (6) and applying the usual boundary-layer approximations: 


ol — mo) OT 
u nD» =" i} (7) 
Ox oy oy oy 7 oy { 


The problem at hand consists in finding solutions to the sys- 
tem of differential equations constituted by Eqs. (1)-(4), to- 
gether with Eq. (7), which satisfy appropriate boundary condi- 
tions. An important simplification is achieved by assuming that 
u, T, and my are functions of a similarity parameter 7 alone, 
where » ~ yx! This assumption is justified as long as the 
boundary conditions permit all relevant quantities to maintain 
constant values along the wall, excepting the normal velocity v. 
The latter may, in general, be permitted to vary inversely as the 
square root of x. When these conditions are satisfied, it is pos- 
sible to reduce the given system to a set of three ordinary dif- 
ferential equations, which is more amenable to integration. The 
transformation is accomplished by means of the following well- 
known substitutions: 

po OY _ OW 


p oy ‘ p Ox 
- (8) 
1 Polto Moltox 


| 
Y=4 f(n) 
2 MoxX \ Po 


where the subscript zero indicates free-stream conditions. 

The following method of solution is singularly suited to the 
present purpose: The set of ordinary differential equations re- 
ferred to above is transformed into a corresponding set of in- 
tegral equations, which are solved by the method of successive 
approximations. It should be pointed out that this procedure 
has been utilized in several reports dealing with boundary-layer 
theory, notably in an interesting paper by H. Schuh.‘ The 
latter has also considered special aspects of the diffusion prob: 
lem, being mainly concerned with the evaporation of liquids in an 
airstream. 

The integral equations expressing conservation of x-momentum, 
conservation of energy, and conservation of foreign gas can be 
written in the general form 


um = [ttm Jo + Kn 


” 
WV, dn + 
0 


Vn dn dn (9 


3) : 
4) 
5) 
| 
3) 
e 
+ 
(eT) 4 (c,T) (« or 4 ou \? 
pu Cp = pe ) 
oy oy \" ay) * “Nay 0 0 
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Fic. 1. Effect of a helium diffusion field at Mach Number 8. 


where [um]o and K,» are constants of integration, and u» repre- 
sents the dependent variables f, (7/7), and mio for m = 1, 2, 
and 3, respectively. In terms of this notation one obtains the 
following expressions for V,, and F,,: 


v= of inl (10) 


Fi = 0 


\ ” 
WV. = ¢3 f, | gaf 


F, = —Prees™ {(1/4) (vo — 1) + (11) 
feu'ToT} 
Vs; = exp leaf dn! 
0 \ (12) 
F, = fnio — — 


where primes indicate differentiation with respect to n, 


41 = g2 = po 'p, gs = kok 

= Cp, Cp, 9s = No 'pomo (13) 

vo =m 
are functions of mo and 7, and Pr denotes the Prandtl Number. 
Comparing this development with the theory of Schuh,‘ one finds 
that Eq. (10) agrees exactly with his results, while Eq. (11) con- 
tains an additional term expressing the effect of variable specific 
heat on the convection of enthalpy. It should be noted also 
that, in the absence of temperature gradients, Eq. (12) loses much 
of its complexity and reduces, in fact, to the integral form of the 
diffusion equation given by Schuh. 

To illustrate the application of the theory outlined above, let 
us consider the following particular problem: A flat plate is im- 
mersed in an airstream at Mach Number 8, and the ratio of plate 
to free-stream temperature is 6. The velocity and temperature 
profiles’ corresponding to this condition are indicated in Fig. 1 by 
broken lines. Now let it be supposed that a foreign gas, which 
we take here to be helium, is introduced into the boundary layer 
near the wall, and that this is done without appreciably altering 


the dynamic and thermal boundary conditions. This require- 
ment is imposed in order to avoid extraneous disturbances that 
would partially mask the primary effect. We stipulate the fol- 
lowing boundary conditions for the concentration profile: nj = | 
atn = 0; mo = Oatn = &. The velocity, temperature, and 
concentration profiles corresponding to the given conditions have 
been calculated with the aid of relevant empirical data® 7 and 
are indicated in Fig. 1 by unbroken lines. It should be noted 
that the direction of heat flow at the surface has been reversed. 
We hope, in a future publication, to present a fuller account of 
the theoretical development, together with a discussion and inter- 
pretation of the numerical results. 
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Experiment on the Decay of Turbulence 
Behind Two Grids 


Hiroshi Tsuji and Francis R. Hama 

Institute of Science and Technology, University of Tokyo, and lowa 
Institute of Hydraulic Research, State University of lowa, 
Respectively 

August 3, 1953 


— EFFECT OF THE UPSTREAM TURBULENCE on the decay of 
turbulence behind a grid was first mentioned by Taylor.! 
He also predicted a decay law based upon the superposition as- 
sumption.? Atsumi’ measured the decay of turbulence behind 
two grids and found essentially no effect of the first grid. The 
writers suspect that the mesh length of the first grid, 14, and 
that of the second grid, M2, are of the same order of magnitude 
and/or that the distance between two grids, x;, is too large in 
Atsumi’s experiment. Several years ago one of the authors 
(Hama‘*) carried out a similar experiment, the mesh length of the 
second grid being much finer than that of the first one. A re- 
markable effect of the upstream grid was observed, especially 
when x;/1/, was less than 50. The turbulence was measured by 
means of the heat-diffusion method,> because a hot-wire ane- 
mometer was not available to the writers at that time, and the 
result was not entirely satisfactory. Recently, the importance 
of this type of experiment has been reemphasized by Lin® 7 and 
Goldstein’ in connection with the validity of the existing decay 
laws. 

The intensity of turbulence behind two grids was again meas- 
ured in the 1.5-m. wind tunnel at the Institute of Science and 
Technology, University of Tokyo, this time by means of the hot- 
wire anemometer. A wooden duct (7.6 m. long and 1 by 1 m. 
square section) was installed to eliminate the free boundary con- 
tamination from the original open jet and to have a longer test 
section. A 0.0003-in. Wollaston wire was welded to the tips of 
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turbulence ahead of the grid. This experiment also confirmed 
ju, Io | / the result of the previous experiment‘ in which Taylor’s super- 

| u2 | / position assumption was found to be invalid. 


gt o> 


the two needles, spaced 5 mm. apart, so as to form an arc, and the 
center part of the arc was etched. The etched part of the wire 
was approximately 1 mm. long. The amplifier had good com- 
pensation characteristics up to approximately 6 ke. for such a 
wire (the time constant was about 1.5 millisec.). The mesh 
length and rod size of the first grid (half-morticed square rods) 
were 5 and 1 cm., and those of the second grid (interfaced circu- 
lar rods) were 1 and 0.2 cm., respectively. The distance between 
the two grids was 50, 100, 150, 200, 250, 350, or 450 em. Under 
this setup the effect of turbulence ahead of the second grid on the 
decay law was expected, since an intense, slowly decaying, large- 
scale turbulence produced by the first grid may survive down- 
stream of the second finer mesh grid. The turbulence level of 
the wind tunnel was approximately 0.2 per cent. 

A summary of the experimental results for the mean velocity of 
10 m. per sec. is shown in Fig. 1, in which x is the distance behind 
the second grid, U is the mean velocity, and u? is the mean- 
square turbulence. For comparison, the decay behind the first 
and second grid alone are also given. Each point plotted is an 
arithmetic mean of at least four measurements from different 
runs. It seems that the turbulence produced by the second grid 
with finer mesh decays rapidly at first. However, the effect of 
the first grid is marked if the distance between the two grids is 
small. As suggested by Lin and Goldstein, the linear decay law 


u? ~ x~! does not hold at the presence of a strong large-scale 


(X-0.06 
0 0.5 1.0 1.5 2.0 25 
Fic. 2. 


in order to check the decay law 
u? ~ [(x — xo)~! + const. ] 


which was originally obtained by Lin® and later endowed with 
firmer reasoning by Goldstein,’ the same data as in Fig. 1 are re- 
plotted in Fig. 2 in a different manner. From Fig. 1, the origin 
Xo is chosen as 0.06 m. Lin’s formula, which is supposed to be 
valid either for the single grid or with a superposed disturbance 
of low frequencies, is confirmed unless the large-scale turbulence 
produced by the first grid is overwhelmingly predominant. The 
discrepancy of the decay curve for the second grid alone may 
have resulted from the relatively high background turbulence. 
The spectrum measurement is now undertaken and will be re- 
ported shortly. 

This investigation was initiated by the suggestion of the second 
author, and, with his close cooperation, ‘the experiment was 
carried out by the first author. The authors would like to thank 
Yoshio Onda for his assistance in conducting the experiment and 
Drs. S. Corrsin, L. S. G. Kovasznay, and C. C. Lin for their help- 
ful suggestions. 
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Some Integrated Properties of Pressure Fields 
About Supersonic Wings 


Z. O. Bleviss 
Aerodynamics Engineer, Douglas Aircraft Company, Inc., 
Santa Monica, Cali 


August 20, 1953 


AND VAN Dyke! have proved some remarkably 
simple theorems concerning the integrated properties of the 
linearized pressure fields on planar supersonic wings, subject to 
certain restrictions on the wing plan forms. It has been found 
that similar properties exist for the pressure fields in certain 
planes about supersonic wings. The results are presented here 
briefly, and one practical application is mentioned. In the follow- 
ing discussion it is to be understood that the results for the inte- 
grated properties of the pressures on the wings have been obtained 
from reference 1. 

Consider a wing with supersonic edges, with its trailing edge 
perpendicular to the flow direction and with an arbitrary dis- 
tribution of local angle of attack a (see Fig. 1). Consider the 
integral of the pressure perturbations p — p.,, in any vertical plane 
(ABCD) which is perpendicular to the wing and parallel to the 
flow direction (note that the vertical plane is not a solid wall); 
the trailing edge (AB) of the vertical plane must be perpendicular 
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to the flow direction and the plane extends forward and vertically 
as far as is necessary to include all pressures forward of AB 
This integral has the dimensions of a force, and the nondimen- 
sional coefficient formed with this force will be denoted by Cp. 
The forward Mach lines from A (AG and A//) and the wing lead- 
ing edge enclose the region of the wing (AGDIH), which pro. 
duces the pressures (dotted shading) in the vertical plane. E de- 
notes the center of pressure (c.p.) of the pressures produced on 
the wing by the angle of attack distribution in region ADH; it 
need not lie in this region. FE is rapidly obtained from the fact 
that it is also the centroid of the angle of attack distribution in 
region ADH. EF is the rearward Mach line (on the wing) from 
E which intersects the vertical plane (or its fore or aft extension) 
at F. Although not indicated, points similar to E and F are as- 
sumed for region AGD. 

Theorem I.—The coetticient Cr has one-half of the two-dimen- 
sional value when based on the area of the region AGDH and the 
average angle of attack in this region—.e. 


Qay 
V 1 V 1 S Seon 


where S = area of AGDH. The coefficient Cy due to any por- 
tion of the region 4GDH is also given by the foregoing statement 
if in that statement the region AGDH is replaced by the portion 
in question. 

Theorem II.—The streamwise position of the c.p. of the pres- 
sures produced in the vertical plane by region ADH is at F. A 
similar statement holds for region AGD. The resultant stream- 
wise position of the c.p. in the vertical plane is obtained by com- 
bining these results using Theorem I. 

The lateral position of the c.p. in the vertical plane is not 
simply related to the wing angle of attack distribution and geom- 
etry. 

The proofs of these theorems follow at once, by superposition, 
from the theorem for the deflected element which is given below. 
A flat-plate wing at zero angle of attack with a deflected element 
producing pressures in the vertical plane (dotted shading) is illus- 
trated in Fig. 2. The trailing edges of the wing and the vertical 
plane are perpendicular to the flow direction. The dimensions of 
the defleeted element are assumed to be infinitesimal (for the 
proper interpretation of such an element see the detailed discus- 
sion in reference 1). The c.p. of the pressures produced on the 
wing by the element is at the centroid (£) of the element. 

Theorem III.—The coefficient Cy has one-half of the two- 
dimensional value (i.e., Cp = — a/0/ M? — 1) when based on the 
area of the deflected element. The streamwise position of the 
c.p. in the vertical plane is at F. 

An extension of Theorem III is immediately obtained from the 
fact that the deflected element may be interpreted as a supersonic 
sink, and, hence, the disturbance field is axisymmetric about the 
streamwise axis through the element. Then the theorem holds 
also for nonvertical planes obtained by rotating the plane 4 BCD 
about this axis. 

The proof of Theorem III was originally obtained by a direct 
integration of the pressures in the vertical plane. However, be- 
cause of the simplicity of the results, a simple superposition proof 
was sought. Such a proof, identical in many respects with that 
given in reference 1, is indicated briefly below. 

Consider any plane above the wing of Fig. 2 obtained by dis- 
placing the wing plane vertically upwards (see Fig. 3). The pres- 
sure field in this plane (dotted shading) is symmetrical about the 
line AD. By the identical superposition procedure used in refer- 

ence | [Figs. 4.1(a) through 4.1(d)], it is easily proved for the dis- 
placed wing plane that (a) the coefficient Cy has the two-dimen- 
sional value —2a0// M? — 1 when based on the deflected area 
and (b) the c.p. is at F, which is the intersection of the Mach line 
EF (from the center of pressure of the deflected elernent and in 
the vertical plane DEF) with the displaced wing plane. For the 
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half plane ABCD the coeflicient Cp has, by symmetry, one-half 
of the two-dimensional value based on the deflected area, the 
streamwise position of the c.p. is still given by F, and the lateral 
position of the c.p. is unknown. Using the axisymmetrical prop- 
erties of the pressure field, the half plane ABCD may be rotated 
into the configuration of Fig. 2 and the theorem is proved. 

One immediate practical application of these theorems is in the 
calculation of interference between fins of cruciform wings. For 
example, if ABCD (Fig. 1) represents a solid wall, the integrated 
properties of the pressures on the side of the wall which is shown 
may be stated as follows, taking proper account of the boundary 
conditions imposed by the wall: The coefficient Cy has the two- 
dimensional value when based on the area of region ADH and the 
average angle of attack in this region and the streamwise position 
of the c.p. is at F. The theorems considered here were, in fact, 
first suggested by numerical calculations of such interference made 
by K. Walker, Jr., who, together with P. A. Lagerstrom, contrib- 
uted further to this work in discussions with the writer 
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Subsonic Oscillatory Aerodynamic Coefficients 
Computed by the Method of Reissner and 
Haskind* 


Gertrude Blancht and Henry E. Fettist 

National Bureau of Standards and Wright Air Development Center, 
Respectively 

August 3, 1953 


7. PROBLEM OF DETERMINING the forces and moments acting 
on an airfoil oscillating in a subsonic flow (when the dis- 
turbances are sufficiently small to permit linearization of the 
equation) can be solved in the two-dimensional case with the aid 
of Mathieu functions. Such a formulation was first given by 
Reissner and Sherman in 1944.' More recently, Haskind? gave 
a more detailed formulation, and Reissner* gave a similar treat- 
ment, using Haskind’s method in part and extending the method 
to include the three-dimensional case. 

Numerical computations using Mathieu functions, but dif- 
fering considerably in detail from the above, have already been 
made by Timman and van de Vooren.‘ Earlier still, Dietze° 
computed tables of lift and moment coefficients based on the 
Possio Integral Equation, in which the known solution for the 
incompressible case was used as an initial approximation. The 
latter method is presumed to be reasonably satisfactory for small 
values of the reduced frequency. Also, independent calculations 
by Turner and Rabinowitz,® using essentially the same method as 
Dietze, gave values that agreed with Dietze’s to within 0.2 per 
cent in magnitude and 0.06° in phase. In addition, one of the 
authors, employing a method that used a polynomial representa- 
tion of the kernal, obtained results which were also reasonably in 
agreement with Dietze’s. The values of Timman and van de 
Vooren,‘ on the other hand, although of the right order of magni- 
tude, departed considerably from Dietze’s, and it was this dis 
agreement that led to the present investigation. 


* This work was conducted at the Institute for Numerical Analysis of the 
National Bureau of Standards and was sponsored by the Aeronautical Re 
search Laboratory, Wright Air Development Center, Wright-Patterson Air 
Force Base, Ohio, 

t Assistant to the Director, Numerical Analysis, Institute for Numerical 


Analysis. 
t Mathematician, Aeronautical Research Laboratory 
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It has been shown by one of the present authors that the formu- 
lations of Haskind? and Reissner* are identical, apart from dif- 
ferences in notation and coordinate systems. Moreover, after 
the Haskind’s paper** was corrected (see reference 8), it was 
found that the method suggested by Haskind is a feasible one for 
numerical computation, especially if further theory of Mathieu 
functions is used to simplify some of the formulas. Since details 
concerning the formulation and method of computation, to- 
gether with a list of errors in references 2 and 3, can be found in 
reference 8, only the final numerical results of the computations 
are included in the present paper. These are contained in Table 
1 which gives the values of the coefficients L,, La, Wi, We for 
Mach Number 0.7. The above coefficients are so defined that 


(as in reference 9) 
total lift = rpb3y?{L,[h/b| + Lala}! 


total moment about forward quarter chord = 
My [h/b] + Malal{ 


where 
h = vertical displacement of forward quarter-chord point 
(positive down ) 
@ = rotation about quarter chord point (positive stalling ) 
p = air density 
b = local semi-chord 
y = forward velocity 
vy = frequency of oscillation 
w = reduced frequency = vb/v 


The values of w/0.02(0.02) (0.10) (0.1) (0.7) | were chosen so as 
to coincide with those of Dietze. Later, two additional values, 
w = 2.914 and 3.643, were taken to determine the behavior of the 
coefficients for larger values of the reduced frequency. The latter 
values were selected so as to cause the parameter u (reference 3) 
to be a convenient quantity. Also included in Table 1 are the 
comparable values computed by Dietze, as well as those of Turner 
and Rabinowitz. Contrary to expectations, the original values 
of Dietz are, in general, closer to the present values than are those 
taken from the later reference. However, both sets of values 
based on Dietze’s method are seen to agree satisfactorily with 
the ones based on the exact solution. Comparison of the moment 
coefficients is somewhat hampered by the fact that, in many cases, 
the original data of Dietze were given to only one significant figure. 
(In fact much closer agreement could have been obtained in many 
cases merely by adding more figures to the original data of Dietze 
with no change in the ones already present. ) 


APPENDIX (I)..-REGARDING THE COMPUTATION OF CONTROL- 
SURFACE COEFFICIENTS 

From the discussion (in reference &) it can be seen that the 
major portion of the labor involved is in the computation of the 
various coefficients in the series determining the pressure distri- 
bution function. Once these are known, it is a comparatively 
simple matter to integrate over only a portion of the chord to ob- 
tain forces and moments on the control surface due to the # and 
a motions. Thus with the information now available, it will be 
possible to obtain values for the coefficients ?,, Pa, 7, Ta which, 
analogously to L and A/ are so defined that 


total force on control surface = wpb v2, Palh/b| + Palal} 
total moment about control surface = Ti [h/b + Ta 


Ordinarily those coefficients would be needed only for values of 
the parameter e from —0.5 to 1.00 (as in reference 9) correspond- 
ing to control-surface-to-chord ratios ranging from 0.75 to zero 

However, because of recently developed relations between the 
above coefficients and those due to control surface motion,'® it is 


** The most serious error in Haskind’s paper was the omission of the nega- 
tive sign in Eq. (2.29). This error cannot be deduced directly from Has 
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Comparative \ 


to 


Th 


= 


: 


x 


Mh: 


Ma: 


by 


w 


Column A* 
Column B 
Column C 


0995 
—0.5162 
2132 
—0.05034 
0.03390 
0.09410 
0.05595 


—43.161 
—18.025 
—9.9905 
—6.4819 
—4.6446 
—0.21993 
—0.13871 


. 2952 
2207 
1577 
. 1057 
.0630 
. 93461 
87472 
83940 
.81069 
7896 
. 73877 
. 10817 
.05409 


— — 


coos 


—6.0649 
—4.2244 
—3.1732 
— 2.4836 
— 2.0009 
—0.89727 
—0.55299 
—0.43442 
—0.40831 
—0.42508 
—0. 45916 
—0.08133 
— 0.08856 


* The values in Column A are correct to four significant fig 


— 129.2 


— 


—36. 


—25 


—19. 
-4. 


TABLE 1 


"alues of the Oscillatory Aerodynamic Coefficients 


Values by method of Reissner and Haskind 
Values by Dietze 
Values by Turner and Rabinowitz 


whenever these functions are nu merically less than unity. 


59994 


.51207 


.14031 


20388 


. 23874 
. 25863 


27025 


. 28830 
30086 
32447 


35917 
40171 


.44706 
35372 


29629 


-918 
242 
. 226 
861 


.935 


5118 
6035 
2173 
4022 
8518 
4388 
33780 
26469 


. 3000 
.2125 
.1556 
1094 
0560 
9305 
. 8739 
8408 
8110 
7792 


.7400 


.0500 
.2188 
.1611 
.4672 
9990 
8995 
5559 
.4338 
. 3990 
4181 


5.000 


175 


.1139 
.4234 


994 
89725 
5561 
4381 
4138 
4258 
4691 
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—17.824 ~17.775 —129.42 —17.7500 —129.25 

—12.515 —12.538 — 59.525 —12.5125 —59.525 
—9.4182 7652 —9.4417 —36.778 —9.3833 —36.75 
—7.3519 779 —7.3781 —25.797 —7.3531 —25.781 

—5.8817 145 —5.8920 —19.410 —5.886 —19.44 
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noted that, if the range of e is extended to include the entire in- 
terval (—1, 1), then coefficients of total lift and moment due 
to control surface motion can be found without the direct con- 
sideration of the corresponding downwash distributions. Thus 
if (as in reference 9) the subscript z denotes vertical translation 
of the control surface and 8 denotes rotation about the control 
surface leading edge, then 


L{e) = Ll, — P,{—e) 
MAe) = My — [1 — (i/w)]Pa(—e) + Pal —e) 
Lae) = [(1/2) — — (i/w)LAe) + T,(—e) — M, 
Male) = Ma — {(1/2) +e — (i/w)] Ma 
(i/w) Mie} + [1 — (¢/w)] Ti(—e) — Tal —e) 


Although it is not possible to compute directly by this method 
the coefficients P., 7., Pg, Tg, the values that can be obtained in 
this way will serve as an indication of the accuracy of those values 
of control surface coefficients now available. 
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Interpretation of Viscous Terms in the 
Turbulent Energy Equation 


Stanley Corrsin 
ee Department, The Johns Hopkins University, Baltimore, 


August 20, 1953 


M*: RECENT PUBLICATIONS using the (correct) turbulent 
kinetic energy equation have misinterpreted the physical 
nature of the viscous term.'~* The equation is due to Reynolds,’ 
who gave both form and interpretation correctly. The only re- 
cent papers in which I have encountered proper interpretation 
without the simplification of local isotropy are two by Rotta,™ !! 
assuming the boundary-layer approximation. 
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The turbulent energy equation for steady-state incompressible 
viscous flow is 


— — -— (ujg*) 
2 ‘ax, ax; 
1 op 
+ (1) 
p OX: 
where 
U; = total velocity 
U; = mean velocity 
u; = U,; — U; (velocity fluctuation) 
= 
vy = kinematic viscosity 
p = density 
p = Static pressure fluctuation 
Xx; = space coordinate 


Physical interpretation of the first four terms is commonly ac- 
cepted as: (a) rate of convection of turbulent energy by mean 
motion; (b) rate of production of turbulent energy from mean 
motion kinetic energy; (c) rate of convection of turbulent energy 
by turbulent motion; and (d) rate of transfer of turbulent energy 
by the work of fluctuating pressure gradients. They need not 
be justified here. 

If we regard Eq. (1) as a transport equation, the viscous term 
can consist of only two possible effects: (e) rate of dissipation of 
turbulent energy to heat (¢) and (f) rate of transport of turbulent 
energy by viscous forces (7)—.e., 


= T — 9 (2) 
The general expression for dissipation rate per unit mass is!* 
Ox, \ Ox; Ox; 
With U; = U0; + uj, it is clear that 
Ou; Ou; Oly, 
Ox; (= =) 


so that Eq. (2) gives the viscous transport as 


ou; Ou; Ou, 
T = vu; 4+ 
OX OX; Ox, \ ON, Ox; 
which can be simplified to 


T =» Vig? + (uae) (5) 


The error committed in references 1 through 8 is that, after the 
mathematically convenient transformation 


= (2 (6) 
i ;= 
4 q Ox; OX, 


it was assumed that on the right side the first term is the full 
transport (because it is expressible as the space derivative of a 
mean quantity?) and that the second term is therefore the full 
dissipation. 

When dissipation occurs in the locally isotropic part of the 


spectrum, 


(2")(2"*) 
Ox; OX, Ox; Ox; 


and the error is negligible. 
Furthermore, when Eq. (1) is integrated across a fully de- 
veloped channel or pipe flow,® the integral of T is zero and the in 


tegral of 


ou, ou k 
= 
Ox; Ox; OX; OX; 


in ¢ is also zero. 


| 
| 
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An alternative interpretation of Eq. (2) follows more closely the 
general discussion in reference.!* 

The viscous term, velocity times viscous force, must be just 
the net rate at which turbulent energy is increased by viscous 
forces. This term must therefore be equal to the difference be- 
tween fofal rate at which viscous forces do work on the volume 
element (1J’) and rate of dissipation to heat—i.e., 


yuu, = W- (2a) 


In other words, W = T. 
I should like to thank Mr. Alan Kistler for his helpful discus- 


sion. 
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The Prediction of Nonlinear Normal Force and 
Pitching Moment for Bodies of Revolution 
Without Boat-Tail 


Howard R. Kelly 
U.S. Naval Ordnance Test Station, Inyokern, China Lake, Calif. 
August 17, 1953 


—— PREDICTION OF NORMAL ¢ORCE AND PITCHING MOMENT for 
bodies of revolution at large angles of attack has never been 
carried out with complete success. One of the most successful 
methods proposed to date was proposed by Allen,'! who assumed 
that the viscous normal force due to the cross flow is independent 
of the axial flow, and may be calculated from the drag coefficient 
Cpe of a cylinder in a fully developed flow normal to its axis. 
Adding this to the slender-body potential theory of Munk,? Allen 
obtains equations that simplify to 


| Ap 
Cup = + (1 — x¢)a? (1) 


for moment coefficient about the base and 
Cy = 2a + nCpc( A p/S)a? (2) 
for normal force coefficient of a body without boat-tail. Allen 


used this Eq. (2) for lift coefficient, but his approximations apply 
better to normal force. In the above equations, 
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FINENESS RATIO 7 CALIBERS 


HEAD LENGTH 2 CALIBERS 7 
| REYNOLDS NUMBER 5.9x 106 

—— COMPLETE THEORY 


—— ALLEN FORMULA 


FOR Cu, 
© DATA FOR Cy 


| 


T 
PITCHING MOMENT COEFFICIENT Cmy 


NORMAL FORCE COEFFICIENT Cn 
\ 


2 Z »4 
! 2 
ar) 5 10 15 20 25 30 
ANGLE OF ATTACK, ow, DEGREES 
Fic. 1. Coefficients of normal force and of pitching moment 


about the base of a tangent-ogive-cylinder at AJ = 0.26. 


a = angle of attack 

n = correction for finite cylinder length 

V = volume of body of revolution 

S = maximum (or base) area of body 

1, = plan form area of body 

1 = length of body in calibers (fineness ratio) 


Xe = axial position of centroid 


The above equations have been found to be inadequate for sub- 
sonic data obtained by the Naval Ordnance Test Station in the 
Consolidated Vultee Aireraft Corporation wind tunnel at San 
Diego. Allen’s theory has been carefully examined, and a new 
method, based on his fundamental equations, has been developed. 
A preliminary version of this method has been published,* and 
further work is in progress to extend the method to drag prediction 
for large angles of attack and to find the upper limit of usefulness 
of the method with respect to angle of attack and Mach Number. 

The principal difference between the present method and that 
of Allen is that the cross flow is here considered not to be fully 
developed but in a state of transient development along the body, 
corresponding to measurements by Schwabe.‘ The results of 
Schwabe for the drag on a circular cylinder started impulsively 
from rest were approximated by a polynomial in odd powers of 
(x tan a)/r, which represents his time parameter. Integration 
along the body shows the increment in normal force coefficient due 
to viscous cross flow is of the form 


ACyw = (2nCp- 'S) (Ay'a® + By'a® + Cy’at +...) (3) 
while a similar equation 
ACy = (2nCp- (A + By'a + Cy’a? +...) (4) 


holds for pitching-moment-coefficient increment. Other modi- 
fications discussed in reference 3 include: 

(a) Showing that Allen’s assumption Cp, = 1.2 is valid only 
for laminar axial boundary layer. If the axial boundary layer is 
turbulent, Cp. = 0.35 gives excellent results. 

(b) Setting » = 1 for bodies without boat-tail. 

(c) Replacing Munk’s potential theory with the second-order 
potential theory of Van Dyke’ for supersonic Mach Numbers. 

(d) Considering the effect of displacement thickness of bound- 
ary layers, especially for turbulent layers in supersonic flow. 

The subsonic equations are found to simplify to 


Cy = 2a + 0.49 (Cp,l?/S)a3 (5) 
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° © DATA FOR Cu, WA = 
Leet © DATA FOR Cy | 6 
° | 
w ° 
as Z 
° 
W 
= 
= = 
° 2 a 6 8 10 
ANGLE OF ATTACK, a, DEGREES 
Fic. 2. Coefficients of normal force and of pitching moment 
about the base of a tangent-ogive-cylinder at AJ = 2.87. 


Cu = (2V/S)a + 0.49( Cp,13/3.S)a3 (6) 


as long as the parameter 2/ tan a does not exceed 4 or 5. Larger 
values of this parameter or extension to supersonic flow will re- 
quire modifications, as mentioned above. 

A typical subsonic case is illustrated in Fig. 1. The boundary 
layer is almost certainly turbulent, so Cp, = 0.35 was used. 
Eqs. (5) and (6) are not adequate for the highest angles of attack 
shown, and curves obtained from the more complete theory are 
also shown in Fig. 1. An example of laminar axial boundary 
layer in supersonic flow is shown in Fig. 2, with data obtained in 
the Naval Ordnance Laboratory tunnel. Eqs. (5) and (6) are 
here modified by using the results of Van Dyke’s theory for the 
linear term. Even better agreement is found if boundary layer 
thickness is considered, but these curves show that the additional 
correction may often be dispensable for laminar boundary layers. 
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Vibration of Rectangular Plates 


J. V. Nagaraja and S. S. Rao 
National Physical Laboratory of India, New Delhi, India 
August 13, 1953 


- teagesantd METHODS as worked out by Rayleigh and Ritz for 
the frequencies and modes of vibration of thin elastic plates 
involve a considerable amount of mathematical computation, 
and the accuracy of the results depends upon the nature of the 
functions set up to represent the plate deflections. With suit- 


able modifications, matrix process applied to the Ritz’ method, 


FORUM 855 


it is possible to arrive at the fundamental mode and frequency of 
thin plates, without much mathematical calculation. 
The plate deflection is represented by 


r 8 


m=1 n=1 


W(x, y) = 
where X and ¥ are functions of x only and y only. These fune- 
tions satisfy the conditions 


and similarly for and where 
derivatives of Y. 

This way of choosing the functions eliminates the difficulties 
involved in solving the simultaneous equations consisting of the 
amplitude coefficients. The plate deflection equation can be 
put in the matrix form as below: 


and X” are first and second 


1 


By substituting this expression in the equation V and 7, we 
minimize it. Here V and T are potential and kinetic energies of 
the vibrating plate and are equal to 


v= SU (= + 
2 oy? 
Ox? Oy? Ox Oy 


phw? 
al W? dx dy (with standard notations) 


The conditions turn out to be 


= w*ph w2 d. 


where 7 and & run through necessary values. 
They take the form 


a b a b 


Giving a particular set of values to the functions it is possible 
to obtain the modes of vibration. The case of a square plate has 
been worked out and the values obtained are compared with those 
by Dana Young. The functions chosen are: 


(1) Clamped-clamped condition: 
fi = x*[(x — a)? — 2a(x — a) + 4a?] 
(2) Clamped-free condition: 
fo = x?/3 (x — a)? 


(3) Free-free condition: 


+ sin? rx/a 


TABLE 1 


Nature of the -——— Value of I Mode—-~ Functions 


Boundary As cal- As per Per cent- Used for 
Conditions culated Dana Young age Error X and Y 
Cantilever plate 3.53 3.494 1.03 fy and fy 
All edges clamped 36.10 35.900 0.30 fo and fe 
Two adjacent edges 8.16 6.958 17.27 fi and fi 
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A Statistical Analysis of Icing-Flight 
Observations 


Emile Z. Gilutin 
Assistant Project Engineer, American Airlines, Inc., Flushing, N.Y. 
August 24, 1953 


SUMMARY 


Icing-flight data were examined statistically from which an equation was 
derived correlating temperature with liquid water content and water droplet 
diameter. Analysis of the errors involved makes feasible the prediction of 
the temperature at which a given set of icing conditions is most likely to 
occur. 


INTRODUCTION 


CING-FLIGHT DATA WERE ANALYZED with the view of obtaining 
I as full an understanding of actually encountered icing condi- 
tions as published data would permit. The ever increasing de- 
mands made by anti-icing systems upon modern high-perform- 
ance aircraft make such an understanding highly desirable. 
The extent to which this is possible in this type of analysis is 
governed by the number of observations, the degree of correla- 
tion between them, and the extent of the errors involved. 


RESULTS 


Three hundred and twenty-five observations of temperature, 
liquid water content, and water droplet from published data! 
yield the following relationship: 


T = — 0.385D — 4.927L + 481.38 (1) 


where 7 is temperature (°R.), D is water droplet diameter 
(microns), and Z is liquid water content (Gm. per cu.m.). The 
degree of correlation, R, among all three variables was found to 
be approximately —0.34, and the standard error of R, computed 
by the relationship 


S.E. of R = (1 — R*)/\/N 


is found to be approximately 0.05. The observed value of the 
degree of correlation being almost seven times its standard error 
seems definitely significant. 

From the frequency-distribution curve for these data, shown in 
Fig. 1, several important conclusions «re possible regarding icing 
encounters. 

(1) Seventy per cent of the time, the encountered temperature 
will not exceed the value obtained in the equation by more than 
approximately +5°R. 

(2) Eighty per cent of the time, the encountered temperature 
will not exceed the calculated value by more than approximately 
+9°R. 

(3) Ninety per cent of the time, the encountered temperature 
will not exceed the calculated value by more than approximately 
+13°R. 

(4) The standard deviation is approximately +10°R. 

(5) In 50 per cent of the encountered icing conditions, the 
temperature will be that actually calculated. 
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DISCUSSION 

It should be emphasized at this point that Eq. (1) is mot an 
algebraic relationship in the sense that the variables may be 
transposed from one side of the equation to the other. The only 
dependent variable is 7, while D and 1 must first be established 
by meteorological and tactical considerations. If the best mean 
D and the best mean L can be established, then this equation will 
provide the best mean 7 at which the given meteorological condi- 
tion will exist. 

An example of the use of this equation follows: Let D = 15 
microns and = 0.8 Gm. per cu.m., then 7 = 472°R. approxi- 
mately. If we wish to provide, conservatively, against 80 per 
cent of the icing encounters (a reasonable assumption) for the 
given D and L, then T = 472°R. — 9°R. = 463°R. or 3°F. 
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On Determining the Additional Apparent Mass 
of a Wing-Body-Vertical-Tail Cross Section 


Richard G. Summers 
Graduate Student, Graduate School of Aeronautical Engineering, 
Cornell University, Ithaca, N.Y. 


August 7, 1953 


E. BRYSON, IN REFERENCE I, developed a method for deter- 
mining the additional apparent mass of the cross section of 
a wing-body-vertical-tail configuration having a circular body, 
symmetrical wings, and unequal vertical-tail spans and trans- 
lating along the nonsymmetrical axis. 
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Using an approach somewhat different from that used by Bry- 
son, one can show that the expression for the additional apparent 
mass (called inertia coefficient 4,; in reference 1) can be obtained 
in a closed algebraic form, rather than the elliptic-integral form 
Bryson found, for the case of a circular body having a straight 
longitudinal axis and zero-thickness wing and tail. 

Consider a right-hand x, y, 2 coordinate system, with the ori- 
gin at the nose of the aircraft, where x is the longitudinal axis, 
positive toward the tail, and where s is positive upward. As- 
sume that the body, as viewed in stationary transverse cutting 
planes (perpendicular to the x axis), is translating and expanding 
or contracting without deforming from its circular shape as the 
aircraft passes through. Let points in this cutting plane be located 
by the complex variable 

w=y+iz (1) 
where y and s are in the same sense as in the x, y, s system, 2 lying 
along the symmetrical axis through the center of the body. Let 
the nonsteady complex potential in this plane be defined as 

Fiw) = @ + wy (2) 

Since the region between the body surface and infinity is free 

of singularities and F(w) is analytic everywhere in this region, we 


can Write 
f F(w) dw = f F(w) dw (3) 
a ce 


where ¢ is the body contour and ¢ is a contour taken in the same 
sense at a great distance from the center of the body. 


If we define 


F(w) = @ + Ws. = complex potential due to expansion or 
contraction of the body 
F'(w) = @' + i’ = complex potential due to cross section at 


rest, in a fluid moving in the —y direc- 
tion at infinity (steady flow) 


then 
F(w) = Fw) + + w (4) 


Substituting Eq. (4) in Eq. (3), we may write 

$ = Im F'(w) dw + f Fw) dw + f w dv | 

ce 
vay (5) 


By virtue of the Cauchy-Goursat theorem and the fact that 


= a constant, on 
= a constant, on ¢ 


Eq. (5) simplifies to 
oa: = Im f F'(w) dw — f dy (6) 
c2 


By Green's theorem, the last integral on the right is simply 


f zdy = —A = cross-sectional area in w plane (7) 
and Eq. (6) becomes, finally, 
f eas = Im f F'(w) dw + A (8) 
C2 
where 
= —m,,/p (9) 


m,, is defined as the additional apparent mass for motion of the 
cross section in the positive y direction. 

The integral on the right side of Eq. (8) can be evaluated by 

Cauchy’s integral formula, provided one can expand the inte- 


grand, for large w, in the form 


FORUM 


F(w)= >> Baw" (10 
a= 
If a convergent series of this form can be found, then 
f F'(w) dw = 27iB,, (11 
cz 
1 is the coefficient of the w~! term in the expansion. 


where B,, 
Finally, Eq. (8) can be written 


A 
Myy = —2rp ea. = ] (12) 


This analysis was applied to the cross section considered by 
Bryson, by mapping the cross section into a flat plate, writing 
the complex potential for this two-dimensional problem, and 
thereby determining the complex potential F’(w) in the w plane. 
Expanding F’(w) for large w, the coefficient of the w~! term is 
found to be 


2 6," 
- — (i +&*) 
16 26.” 
+ + (1 + (13) 
where 
5 a 5 a 5 a 
and 


b = h(x) = span of wing 


a = a(x) = radius of body measured — radially 
4, = (x) = height of top vertical tail from the center of 
ty = fo(x) = height of bottom vertical the body 
tail 
For the circular body considered, 
A/2x = a?/2 (15) 


and from Eq. (12) the desired formula becomes 


6,2 2 (] 
= — (1 + 5,7)? — 95.2 | + §.4) 
— - — X 
04 
25.2 + &*) + 52! + 62 
| 5,2 
yi + 6) + + 634) (16) 
where m,,, = m,,(X). 


Since B,, is purely real, it can be shown that 


F oay = Re f F'(w) dw = 0 (17) 


which verifies Bryson’s statement that 4, = 0 for aircraft of the 
type considered, having one plane of symmetry. Eq. (16) has 
been applied to several special cases for which results are avail- 
able, and complete agreement has been found. 
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Pressure Rise Required for Separation in 
Interaction Between Turbulent Boundary 
Layer and Shock Wave 


Robert D. Tyler and Ascher H. Shapiro 

Assistant Professor of Mechanical Engineering (Visiting, on Leave 
from Imperial College of Science and Technology, London, 
England), and Professor of Mechanical Engineering, Massachusetts 
Institute of Technology, Cambridge, Mass., Respectively 


August 13, 1953 


INTRODUCTION 


— EXPERIMENTS! ° have shown that shocks of suffi- 
cient strength produce separation of boundary layers. In 
such cases the shock that is directly associated with separation is 
not that which would be present in the absence of a boundary 
layer but is rather that which is induced by the separation itself 
(Fig. 1). Indeed, in the zones marked ‘‘region studied”’ of Fig. 
1, there is a feedback between the shock and the boundary 
layer, inasmuch as the pressure rise for separation is produced by 
the shock, and the shock is, in turn, generated because of the flow 
deflection associated with separation. 


Mope_ EMPLOYED FOR ANALYSIS 


When the boundary layer is turbulent, it is observed! that the 
pressure rise occurs in an extremely short distance, of the order of 
a few boundary-layer thicknesses. This suggests an approximate 
but simple discontinuity analysis for determining the pressure 
rise for separation, based on the assumptions that in the zone of 
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sudden pressure rise the skin-friction force and the mass flow 
brought into the boundary layer from the free stream are both 
negligible. 


ANALYSIS 


Under these conditions the edge of the boundary layer is a 
stream line. Considering the control volume ABCD of Fig. 2, 
the continuity law and momentum theorem may be written as 


Continuity: pu’ dy’ = dy” 
0 0 
Momentum: 
— = dy” — pu’? dy 2) 
2 0 0 


where the pressure force along AB has been linearly approxi- 
mated by (pi + p2) (62 — 6:)/2, and the pressures along 1D and 
BC have been assumed uniform and equal to p, and ps, respec- 
tively. 

Introducing the customary definitions of displacement and 
momentum thickness, 


1 p’ n’ y’ 
= 1 — — 
0 pry 


0 1 p’ u’ 
= = 1 -id = Cae: 
6 1 0 Pl ity 6; 


Eqs. (1) and (2) may be combined and rearranged algebraically 
to give 


2kM 2 
1+ 


pits | 1 — 
p2lte 1- (5*/5)2 


(4) 


Since boundary layers cannot support large sudden pressure rises, the shock induced by separation must be comparatively weak, and 
therefore the entropy change in the free stream between sections AD and BC may be ignored with little error. Considering the external 


flow as isentropic, we may write pi* = p2*, pi* = p2*, and m* = 
flow. With these and the formula 


Mpi/pi* = 


Eq. (4) may be written in the form 


CHOICE OF VELOCITY PROFILES 
The values of 6*/6 and 6/6 were evaiuated with the assumption 
of constant stagnation temperature in the boundary layer. 
At section AD the velocity profile was taken as a power-law 
profile corresponding to turbulent flow with zero pressure gra- 
dient. 


Turbulent: = with 2 = 5, 7,9, 11 (6) 
At section BC the profile was taken to be on the verge of sep- 
aration, with 0p/Ox = 0. Expressing «”/u2 as a power series in 
y”/d2 and noting that the differential equations of the boundary 
layer require the first, second, and third derivatives to be zero 
at the wall and the first derivative to be zero at y"/d2 = 1, we may 
formulate such a profile as 


u.*, where the asterisk refers to sonic conditions in the isentropic 


‘uy ** p,* 


= (8*/8); — (6/8): — (8*/8)2 — (0/8): 
I 1 — (6*/6), 1 — (5*/8)2 


Separating: — =! = 7) 


With the foregoing assumptions, 5*/6 and 6/6 at sections 1 and 
2 may be evaluated from Eq. (3) in terms of 1M, or My. Accord 
ingly, all the quantities in Eq. (5) are functions only of .V, or 
M,. Consequently Wo, the free-stream Mach Number at the 
separation point, may be determined for any value of .,, the 
free-stream Mach Number just ahead of the abrupt pressure 
rise. Knowledge of V2 and JM, leads to the corresponding value 
of 


P = (po — pi)/(1/2) pits? 


the dimensionless pressure rise supportable by the boundary 
layer in an abrupt separation. 
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RESULTS 

The results are presented in Fig. 3, where 7? is plotted against 
M, for different values of n. 

Employing values of » of 7, %, and 11 for length Reynolds 
Numbers of 16°, 10°, and 107, respectively, P is plotted versus 
Re, for fixed values of WV, in Fig. 4. Experimental values® ® 8 
are shown for comparison in Fig. 4. It is gratifying to note that 
the analytical results give correctly the order of magnitude of the 
observed results, thus indicating that the analysis incorporates 
the principal physical phenomena of the process. 
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Measurements® indicate that P decreases slightly with 1/, in 
the supersonic range. The present analysis yields the same re- 
sult. 

The theoretical prediction that ? should increase with increas- 
ing Reynolds Number is opposite to the experimental trend 
shown in Fig. 4. Thus it is evident that the terms omitted from 
Eqs. (1) and (2) (skin friction; mass and momentum introduced 
into the boundary layer) are of some importance. 

Examination of the complete momentum equation of the 
boundary layer shows that the momentum flux brought into the 
boundary layer from the free stream exceeds the skin-friction 
drag. Accordingly, the present analysis should give too low a 
value for P?, inasmuch as the incoming momentum augments the 
pressure rise supportable by the boundary layer. Both "the 
incoming momentum flux and skin friction will vary approxi- 
mately as dé*/dx—i.e., approximately as Re;~'/"~? for a tur- 
bulent layer. This reasoning shows that the effects omitted from 
the analysis would tend to make P? diminish as the Reynolds 
Number increases. It is planned to extend the present analysis 
to take account of wall friction and the entrained free-stream 


momentum. 
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Effects of Wing Twist on the Response of an Airplane Encountering 
a Sharp-Edged Gust 


(Concluded from page 834) 


To estimate the magnitude of this parameter for a given 
airplane, we proceed as follows: 

If we define m, as the slope of the lift coefficient curve 
for the airplane with an elastic wing and m that of the 
airplane with a rigid wing so that the following relations 
exist: 

Cr = ma = m(a — Aa) = m(a +6) (Cl) 

[since, from the assumption Eq. (1), if the elastic center 

lies ahead the aerodynamic center, we should have Aa = 
—6], we must have 

Mz = —bL = —(1/2)pU?Sma (C2) 

= —(1/2)pU*Sm(a + 6) (C3) 


Using the relations, Eq. (28) and Eq. (1), we reduce 
Eqs. (C2) and (C3) to the following forms, respectively: 


1 
— MU? ( ) gts) = pl? Sma (C4) 


2\2 1 
( ) 2"(s) = pU?Sma + 
2 
1 
5 PU*Sm Rz"(s) (C5) 


Eliminating s”(s) from Eqs. (C4) and (C5), we have 


Denoting the following ratio by 7, 
r= m,/m (C7) 
and using the symbol \ defined in Eq. (20), we obtain 


cR = 4\[(1 — 7)/r] (C8) 


{It can be noted that if the elastic center lies behind the 
aerodynamic center, as in most straight wing aircraft, 


Aa = +0andcR = —4A{(1 — r)/r}}. Hence, we con- 
clude that cR is a dimensionless parameter that can 
be expressed in terms of other two dimensionless 
parameters } and ry. Given an airplane, the value of 
cR may be obtained from Eq. (CS) if the parameters 
and r for that airplane can be reasonably estimated. 
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Since its incorporation in October, 1932, as a scientific membership society to advance the art and science 
of aeronautics, the Institute of the Aeronautical Sciences, in over twenty years of service to aviation, has 
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become the representative technical society for the aviation industry and profession. It is the only organiza- 
tion in the United States which brings together all of the sciences and branches of engineering and tech- 


nology which are applied to aeronautics. 


Publications 


Papers presented at meetings or submitted to the Edi- 
torial Board are published in the JourNaL or THE AERO- 
NAUTICAL ScrENcEs or the AERONAUTICAL ENGINEERING 
Review upon recommendation by the Editorial Com- 
mittee. The JourNnat prints in each issue full-length 
scientific papers on new research and developments 
in various fields applied to aviation. The Review keeps 
members and subscribers up to date on aeronautical news 
and literature through reviews of new books and periodical 
articles, government publications and trade literature and 
publishes papers on applied engineering and aircraft pro- 
duction. 

The AgronauticaL ENGINEERING CataLoc, published 
yearly, is a guide to sources and specifications of materials, 
parts, and accessories used in the design and production of 
aircraft, aircraft engines, and parts. 


Research Facilities 


Through the Aeronautical Archives of the Institute, 
members have for their use the most complete aeronautical 
library and reference research facilities maintained by any 
technical society. At the New York building of the Insti- 
tute, The W. A. M. Burden Library of books, reports and 
periodicals, the Aeronautical Index of subject files, bio- 
graphic files and bibliography, and the historical and art 
collections in the Archives are available for study. The 
Pacific Aeronautical Library at 7660 Beverly Boulevard, 
Los Angeles, maintains a reference collection and reading 
room for members and loans books to company libraries 
in the area. Through the Paul Kollsman Lending Li- 
brary, aeronautical books are lent by mail without charge 
to members anywhere in the continental United States. 
In addition to a complete collection of standard aero- 
nautical reference works, the Kollsman Lending Library 
has available for loan new aviation books as they are re- 
viewed in the AERONAUTICAL ENGINEERING Review. 


Membership 


In addition to membership in professional societies 
representing their particular fields, specialists applying 
their knowledge and experience to aeronautics find mem- 
bership in the Institute of great value in the contacts made 
with others engaged in aeronautical work and in the broad 
coverage of all aeronautical problems furnished through its 
meetings and publications. All applications for member- 
ship are carefully reviewed by the Membership Committee, 
which recommends the grade of membership for which the 


INSTITUTE OF THE AERONAUTICAL SCIENCES, INC. 


2 EAST 64TH STREET 


NEW YORK 2], N.Y. 


applicant is found eligible. There are advanced grades of 
membership—Fellow and Honorary Fellow—to which 
members may be elected as their experience warrants. 
Thus, a member receives due recognition—more valuable 
because it is given by his colleagues—for his contributions 
to the progress of aeronautics. New applicants may be 
elected by the Membership Committee to the grade of 
Associate Fellow if they have been engaged for not less 
than ten years in the practice of the aeronautical sciences 
and shall have been in responsible charge of important 
scientific or engineering work, or shall have otherwise 
made outstanding contributions to the aeronautical sci- 
ences; to the grade of MEMBER if they have been en- 
gaged in professional work for at least eight years and have 
acquired a recognized standing in engineering, design, or 
other special work applied to aeronautics; to the grade of 
Associate Member if they have acquired a recognized 
standing in an administrative capacity in the aviation in- 
dustry; to the grade of Technical Member if they are recent 
engineering school graduates or are engaged in technical 
aeronautical work. Student Members are admitted by 
application through Student Branches organized at their 
schools. Aeronautical companies, as well as individuals, 
are affiliated with, and participate in the support of, the 
Institute through Corporate Membership. 


Dues 


An entrance fee of $10 is required of all new members 
who apply for admission to a grade of membership to be 
specified by the Membership Committee (Corporate Mem- 
bers, Foreign Members, and Student Members being ex- 
cepted). 

Annual membership dues include a subscription to the 
AERONAUTICAL ENGINEERING Review. Dues for the vari- 
ous grades of membership are as follows: Fellows— 
$16.50; Associate Fellows—$16.50; MEMBERS—$13.50; 
Associate Members—$13.50; Technical Members—(26 
years of age and over) $13.50, (under 26) $6.50. 

The regular subscription price for the JouRNAL OF THE 
AERONAUTICAL SciENcEs and the AERONAUTICAL ENGINEER- 
1NG Review is $12 and $3.00 per year, respectively. Mem- 
bers may subscribe to the JourNat at $4.00 per year. 

Those who are engaged in aeronautical work or who 
have an interest in any technical phase of aviation will find 
the services of the Institute to be a necessary aid in keeping 
abreast of developments in these times of such rapid prog- 
ress. An application form and further information about 
membership can be obtained from an Institute member or 
by writing to the Secretary. 
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SUGGESTIONS for CONTRIBUTORS 
to the Publications of the 
INSTITUTE of the AERONAUTICAL SCIENCES 


The Institute of the Aeronautical Sciences invites both members and nonmembers from any 
country to submit papers for publication in the JOURNAL OF THE AERONAUTICAL SCIENCES and the 
The-Institute, following the practice of other societies, does 


The following directions for the preparation of papers, if followed by authors, will save corre- 
spondence, avoid the return of papers for changes, minimize the work of preparation for the printer, 
and save the expense due to the charges made for “‘author’s corrections.” 


Manuscrirts: The original typewritten copy of the paper is 
| desired, double or triple spaced on one side of white paper sheets, 
consecutively numbered. There should be wide margins to 
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